AMERICAN

METEOROLOGICAL
SOCIETY

Journal of Climate

EARLY ONLINE RELEASE

This is a preliminary PDF of the author-produced

manuscript that has been peer-reviewed and
accepted for publication. Since it is being posted
so soon after acceptance, it has not yet been
copyedited, formatted, or processed by AMS
Publications. This preliminary version of the
manuscript may be downloaded, distributed, and
cited, but please be aware that there will be visual
differences and possibly some content differences
between this version and the final published version.

The DOI for this manuscript is doi: 10.1175/JCLI-D-16-0225.1

The final published version of this manuscript will replace the

preliminary version at the above DOI once it is available.
If you would like to cite this EOR in a separate work, please use the following full
citation:

Xue, P., J. Pal, X. Ye, J. Lenters, C. Huang, and P. Chu, 2016: Improving the
Simulation of Large Lakes in Regional Climate Modeling: Two-way Lake-
atmosphere Coupling with a 3-D Hydrodynamic Model of the Great Lakes. J.
Climate. doi:10.1175/JCLI-D-16-0225.1, in press.

© 2016 American Meteorological Society

Manuscript (non-LaTeX) Click here to download Manuscript (non-LaTeX) TC-
GLARM_Revision_201609_submit_version.docx

1 Improving the Simulation of Large Lakes in Regional Climate Modeling:

2 Two-way Lake-atmosphere Coupling with a 3-D Hydrodynamic Model of the Great Lakes

4 Pengfei Xue1*, Jeremy S. Pal2, Xinyu Ye1, John D. Lenters3, Chenfu Huang1, Philip Y. Chu4

1
6 Great Lakes Research Center, and Department of Civil and Environmental Engineering,

7 Michigan Technological University, Houghton, Michigan

2
8 Department of Civil Engineering and Environmental Science, Seaver College of Science and

9 Engineering, Loyola Marymount University, Los Angeles, California

3
10 LimnoTech, Ann Arbor, Michigan
4
11 NOAA - Great Lakes Environmental Research Laboratory, Ann Arbor, Michigan

12

*
13 Corresponding Author: Pengfei Xue ([email protected])

14

15 September 9, 2016

16

17

18

19

20

21

22

23

24
25 Abstract

26 Accurate representations of lake-ice-atmosphere interactions in regional climate

27 modeling remain one of the most critical and unresolved issues for understanding large-lake

28 ecosystems and their watersheds. To date, the representation of the Great Lakes two-way

29 interactions in regional climate models is achieved with 1-D lake models applied at the

30 atmospheric model lake grid points distributed spatially across a 2-D domain. While some

31 progress has been made in refining 1-D lake model processes, such models are fundamentally

32 incapable of realistically resolving a number of physical processes in the Great Lakes. In this

33 study we develop a two-way coupled 3-D climate-lake-ice modeling system (named TC-3D-

34 GLARM) aimed at improving the simulation of large lakes in regional climate models and

35 accurately resolving the hydroclimatic interactions. Model results are compared to a wide variety

36 of observational data and demonstrate the unique skill of the coupled 3-D modeling system in

37 reproducing trends and variability in the Great Lakes regional climate, as well as in capturing the

38 physical characteristics of the Great Lakes by fully resolving the lake hydrodynamics.

39 Simulations of the climatology and spatiotemporal variability of lake thermal structure and ice

40 are significantly improved over previous coupled, 1-D simulations. At seasonal and annual time

41 scales, differences in model results are primarily observed for variables that are directly affected

42 by lake surface temperature (e.g., evaporation, precipitation, and sensible heat flux) while no

43 significant differences are found in other atmospheric variables (e.g., solar radiation, cloud

44 cover). Underlying physical mechanisms for the simulation improvements using TC-3D-

45 GLARM are also discussed.

46

47

48

2
49 1. Introduction

50 Consisting of five large lakes, the Great Lakes of North America (also known as the

51 Laurentian Great Lakes) form the largest surficial area of fresh water on Earth. The total surface

52 area of the Great Lakes covers approximately 245,000 km2, spanning approximately 16 degrees

53 of longitude and 7.5 degrees of latitude (Figure 1). The Lakes contains 23,000 km3 of freshwater,

54 holding one-fifth of the world's surface fresh water by volume. Lake depth varies significantly

55 from a few meters in the coastal regions to a few hundred meters in the deep basins, with the

56 deepest location being roughly 400 m in southeastern Lake Superior (Table 1). Due to their sea-

57 like characteristics (distant horizons, great depths, steep bathymetric gradients, rolling waves,

58 sustained winds, strong Coriolis-influenced currents, and large thermal variability), the Great

59 Lakes have long been referred to as “inland seas,” playing a critical role in the variability of the

60 hydroclimatic system throughout the Great Lakes region [Williamson, 1854] .

61 The Great Lakes are not only sensitive to climate change, which is likely to have

62 contributed to recent large fluctuations in lake level, thermal structure, and ice coverage, but they

63 are also a significant regional climate driver due to their large volume, thermal inertia, and

64 surface area [Wang et al., 2012; Clites et al., 2014a; Van Cleave et al., 2014; Gronewold et al.,

65 2015] . Each of these factors contribute to strong linkages between air and lake temperature,

66 evaporation, precipitation, and ice coverage in the coupled regional lake-atmosphere system

67 [Blanken et al., 2011; Lenters et al., 2013] . Lake surface temperature (LST) and ice coverage

68 have a strong impact on regional climate. They also serve as physical indicators of climate

69 change through direct interactions with surface winds and atmospheric heat and moisture fluxes,

70 which in turn modify the lakes’ thermal structure, water level, and circulation, making the

71 system particular sensitive to climate change. For example, warmer winters, loss of ice coverage,

3
72 and earlier summer stratification of Lake Superior have caused summer LST to warm faster than

73 summer air temperatures in recent decades [Austin and Colman, 2007], while the reverse has

74 been observed in winter, with air temperature warming faster than winter LST [Lenters 2004].

75 Rapid summer LST warming has also been seen in other deep lakes around the world [Huang et

76 al. 2012, O'Reilly et al., 2015, Zhong et al. 2016] . In addition, large interannual variations in

77 winter air temperature and ice coverage on the Great Lakes have been observed recently within

78 only a three-year time span, often referred to as the “big heat” of 2012 and “big chill” of 2014,

79 leading to significant impacts on water temperature, evaporation, and lake levels in the

80 subsequent ice-free seasons [Wang et al., 2012; Lenters et al., 2013; Clites et al., 2014a;

81 Gronewold et al., 2015] . This has brought renewed attention to the future impacts of large trends

82 and strong climatic variability on the Great Lakes. Nevertheless, the complexities of the thermal

83 structure, ice distribution, and circulation patterns in these large “inland seas” make the

84 prediction of such impacts extremely challenging.

85 Climate models are the primary tool used to assess climate change and associated impacts

86 [IPCC, 2013] . Despite the significant effects of the Great Lakes on the regional climate, only

87 four of the Coupled Model Intercomparison Project Phase 5 (CMIP5) Atmosphere Ocean

88 General Circulation Models (AOGCMs) used in the latest IPCC report [IPCC, 2013] provide

89 even crude representations of the lakes (less than 20 grid points for the Great Lakes); the

90 remaining AOGCMs treat the lakes as land points. Regional climate models (RCMs) aim to

91 enhance regional detail in response to regional-scale forcing through a more realistic

92 representation of physics and dynamics by including finer-scale topography, vegetation, and

93 land/water coverage [Feser et al., 2011; Giorgi, 2006] . The finer grid resolution (tens of

94 kilometers) generally improves simulations of regional climate and provides more detailed

4
95 characteristics of temperature, wind, moisture, and precipitation in comparison to global

96 reanalysis datasets and AOGCM simulations. To date, the most common coupled representations

97 of the Great Lakes region are performed using RCMs coupled with one-dimensional (1-D) lake

98 models at RCM lake-grids (i.e., 1-D lake models are distributed spatially across a two-

99 dimensional (2-D) domain to form a three-dimensional (3-D) representation of a lake, hereafter:

100 RCM coupled with 1-D lake models) [Hostetler et al., 1993; Lofgren, 2004; Subin et al., 2012;

101 Notaro et al., 2013] . Although considerable progress has been made in refining the 1-D deep

102 lake model processes, primarily through an improved characterization of vertical mixing and

103 eddy diffusivity [Subin et al., 2012; Bennington et al., 2014; Lofgren, 2014] , Bennington et al.,

104 [2014] point out that “Although the Hostetler lake model within RegCM4 is now able to capture

105 a reasonable vertical structure of temperature, circulation and dynamics must be accounted for in

106 large, deep lakes.” Other studies applying RCMs with 1-D lake models to the Great Lakes region

107 also recognize the need for 3-D lake models to accurately represent the physical characteristics

108 of the Great Lakes and properly resolve horizontal and vertical mixing processes to reduce biases

109 in LST, ice coverage, and thermal stratification in regional climate simulations [e.g., Lofgren,

110 2004; Gula and Peltier, 2012; Notaro et al., 2013] .

111 Studies in other regions also report that traditional modeling approaches that neglect lake

112 hydrodynamics and consider only thermodynamics alone are not satisfactory. For example, Song

113 et al., [2004] found that not including the wind-driven transport of heat from warm regions of

114 Lake Victoria to cooler regions within the lake results in a degraded simulation of the climate

115 downstream, not only over the rest of the lake, but also the surrounding land regions. The

116 resulting asymmetric LST distribution modifies the overlying wind circulation, cloud cover, and

5
117 rainfall. Only a 3-D lake model coupled with an RCM resolves this secondary feature through

118 explicit lake–atmosphere interactions.

119 All of these research concerns point to the urgency of properly resolving large-lake

120 hydrodynamics and interactions with the overlying atmosphere for improving regional climate

121 modeling of the Great Lakes. The most suitable approach is to represent the system through two-

122 way coupling of a regional climate model with a 3-D hydrodynamic model. This is the current

123 direction of next-generation model development and the most likely method for obtaining

124 reliable projections of the future impacts of climatic trends and interannual variability on the

125 Great Lakes system from a regional modeling perspective. Thus, the purpose of this paper is to

126 introduce such a two-way coupled modeling system for the Great Lakes, including

127 demonstration of the reliable skill of the model in simulating large-lake hydrodynamics and

128 regional climate dynamics over the Great Lakes region. The seasonal and interannual variability

129 of the regional climate, lake circulation, thermal structure, and ice cover of each lake are

130 examined, along with estimates of surface heat and moisture fluxes.

131 The remaining parts of this paper are organized as follows. In sections 2 and 3, we

132 describe the configuration of the coupled lake-atmosphere model, the design of the numerical

133 simulations, and the model validation data used in this study. In section 4 we discuss the

134 modeling results in comparison with a wide array of in situ and satellite datasets. The results of

135 two extreme events (the 2012 “big heat” and 2014 “big chill”) are examined in section 5,

136 followed by discussion and conclusions in section 6 and 7.

137

138 2. Models and data

6
139 To the best of our knowledge, this is the first documentation of a two-way coupled, 3-D

140 regional climate modeling system for the Great Lakes, including a 3-D hydrodynamic model of

141 lake circulation, thermal structure, and ice dynamics. For brevity, we hereafter refer to the new

142 modeling system as the “Two-way Coupled 3-D Great Lakes-Atmosphere Regional Model,” or

143 TC-3D-GLARM to distinguish it from previous models that use 1-D, 2-D, and/or uncoupled (or

144 one-way coupled) numerical modeling approaches.

145 Regional climate model. The latest (4th) version of the International Centre for

146 Theoretical Physics (ICTP) Regional Climate Model (RegCM), known as RegCM4, is used to

147 simulate land and atmospheric processes [Giorgi et al., 2012] . RegCM4 is a 3-D, hydrostatic,

148 compressible, primitive equation, σ-coordinate regional climate model based on the hydrostatic

149 version of the fifth-generation Pennsylvania State University–National Center for Atmospheric

150 Research (PSU–NCAR) Mesoscale Model (MM5) [Grell et al., 1994; Pal et al., 2000] . The

151 model physics are similar to those of RegCM3 and are described in detail in Pal et al., [2007],

152 Steiner et al.,[2009], and Wang et al. [2015] . It adopts the Community Climate Model Version 3

153 (CCM3)-based package for atmospheric radiative transfer computations [Kiehl et al., 1996] ,

154 and nonlocal formulation of the planetary boundary layer from Holtslag et al., [1990] . The

155 large-scale cloud and precipitation schemes are resolved by the explicit cloud/precipitation

156 scheme of Pal et al. [2000] , and the unresolvable precipitation processes (cumulus convection)

157 are represented using the Grell parameterization [Grell, 1993; Grell et al., 1994] . The model

158 version used in this study is based on Wang et al. [2015] with carbon-nitrogen dynamics (CN)

159 and dynamic vegetation (DV) components turned off, and surface physics calculations are

160 performed using the Community Land Model 4 (CLM4) to represent soil-vegetation

7
161 hydrological processes of the land surface. Surface fluxes over water (e.g., oceans and lakes) are

162 handled by Zeng’s bulk aerodynamic ocean flux parameterization scheme [Zeng et al., 1998] .

163 RegCM4, by default, can be coupled to the 1-D, energy-balanced, diffusion-convection

164 lake model described by Hostetler et al., [1993] , as well as the one-layer ice model from

165 Patterson and Hamblin [1988] at each model lake gridcell. In this configuration, the 1-D models

166 are distributed spatially across a 2-D domain, and the energy is transferred between layers by

167 eddy and molecular diffusion and by vertical convective mixing. This coupled RegCM4 1-D

168 lake-atmosphere modeling system has been applied to the Great Lakes in a number of previous

169 studies [Notaro et al., 2013; Bennington et al., 2014; Notaro et al., 2015] . In particular,

170 Bennington et al., [2014] demonstrated in great detail the performance of the coupled RegCM4

171 1-D lake simulation for the Great Lakes, including challenges and failures with the default

172 model, improvements in the modified version, and the overall, fundamental limitations of such a

173 model.

174 Our RCM modeling domain covers not only the Great Lakes region but also the majority

175 of North America. There are 350 x 360 horizontal grid points at 18-km grid spacing (Figure 1)

176 and 18 vertical sigma layers. It is worth noting that most previous RCM studies for the Great

177 Lakes have been isolated to the immediate surrounding region, which can constrain the system to

178 the large-scale, driving AOGCM dynamics. These studies have also been performed at a coarser

179 resolution (e.g., 25 to 60 km), which can leave many finer-scale processes unresolved. The

180 lateral atmospheric boundary conditions for our regional climate modeling system are provided

181 by ERA-interim climate reanalysis data from the European Centre for Medium-Range Weather

182 Forecasts (ECMWF). Lateral boundary conditions include 6-hourly surface pressure and wind

183 components, air temperature, and mixing ratio at all vertical model levels.

8
184 Hydrodynamic model. The Great Lakes hydrodynamic model used in this study is based

185 on the Finite Volume Community Ocean Model (FVCOM) [Chen et al., 2006] , a free-surface,

186 primitive-equation hydrodynamic model that solves the momentum, continuity, temperature,

187 salinity (often not used for the Great Lakes), and density equations and is closed physically and

188 mathematically using turbulence closure submodels. FVCOM is solved numerically using the

189 finite volume method in the integral form of the primitive equations over an unstructured

190 triangular-grid mesh and vertical sigma layers. This approach combines the best features of an

191 unstructured grid for ideal geometric fitting and the flexibility of local mesh refinement (similar

192 to finite-element methods), as well as numerical efficiency and code-simplicity (similar to finite

193 difference methods). FVCOM has been widely implemented in coastal ocean applications [Xue

194 et al., 2011; Xue et al., 2012] , as well as the Great Lakes [Shore, 2009; Anderson and Schwab,

195 2013; Xue et al., 2015]. Other hydrodynamic models have also been used for the Great Lakes

196 (e.g., [Schwab and Bedford, 1994; Beletsky et al., 2006; Wang et al., 2010; Huang et al., 2010;

197 Fujisaki et al., 2012] ), but all of these models were developed in one-way or uncoupled modes.

198 The horizontal resolution of the model triangular grids varies from ~1-2 km near the

199 coast to ~2-4 km in the offshore regions of the lake (Figure 2). The model is configured with 40

200 sigma layers to provide a vertical resolution of < 1 m for nearshore waters and ~2-5 m in most of

201 the offshore regions of the lakes. The Mellor-Yamada level 2.5 (MY25) turbulence closure

202 model [Mellor and Yamada, 1982] is used for simulating vertical mixing processes, which

203 includes a set of prognostic equations for turbulent kinetic energy and a length-scale-related

204 parameter to calculate eddy viscosities and vertical diffusivities. The horizontal diffusivity is

205 calculated using the Smagorinsky numerical formulation [Smagorinsky, 1963], determined by the

206 horizontal velocity shear as well as the model grid resolution. We note here that although the five

9
207 Great Lakes are geographically distinct, Lakes Michigan and Huron are hydraulically connected,

208 so the system can often be treated as four enclosed basins when hydrodynamics and regional

209 climate are of interest.

210 Ice model. To simulate ice-water interaction processes, we adapted the Los Alamos

211 Community Ice Code (CICE) into an unstructured-grid, finite-volume version within the

212 FVCOM framework [Gao et al., 2011] . CICE uses four layers of ice, one layer of snow, and is

213 governed by energy-conserving thermodynamics [Maykut and Untersteiner, 1971; Bitz

214 and Lipscomb, 1999] , elastic-viscous-plastic ice momentum equations [Lipscomb and Hunke,

215 2004] , and energy-based ridging schemes of [Thorndike et al., 1975; Lipscomb et al., 2007] ,

216 with ice strength parameterizations given by [Rothrock, 1975] . Ice cover in the Great Lakes

217 typically lasts 3-5 months every year and is confined within the closed basins. Thus, we

218 developed a simplified and much more computationally efficient version of CICE as an energy-

219 conserving thermodynamic model of ice. Its governing equation is well described by Bitz

220 and Lipscomb [1999] , and we adapted the ice fraction calculation at each individual model cell

221 by calculating the evolution of the ice thickness distribution in time and space through

222 thermodynamic growth and melting.

223 The simplified ice model is very computationally efficient, allowing us to conduct

224 coupled atmosphere-hydrodynamic-ice simulations that span more than a decade in length. Even

225 so, it is still a more sophisticated model than the default ice model in RegCM4, since it allows

226 for ice fraction calculation at each model cell based on the ice thickness distribution function

227 (ITD). More importantly, the model is employed in the FVCOM model grid at a resolution of 1-2

228 km, which is much higher resolution than the ice model employed in RegCM4 (18-km grid

229 spacing).

10
230 Data for Model Validation. To validate the atmospheric components of the TC-3D-

231 GLARM modeling system, we focus on the climatic trends, spatial patterns, and interannual

232 variability in surface air temperature and precipitation, both at the synoptic scale of North

233 America and the regional scale of the Great Lakes. Model results are evaluated against station-

234 based observational datasets over North America, namely the global land-station-based half-

235 degree datasets produced by the Climate Research Unit (CRU, version 3.0) for the period 1982-

236 2013 [Harris et al., 2014] . Over the Great Lakes, we evaluate the model-simulated precipitation

237 and evaporation fields against estimates from the Great Lakes Hydro-Climate dashboard

238 (GLHCD), which are products based on a combination of models and observations [Gronewold

239 et al., 2013; Clites et al., 2014b; Hunter et al., 2015] .

240 For comparisons of LSTs, we evaluate the model results against the Great Lakes Surface

241 Environmental Analysis (GLSEA2), which is available from

242 https://coastwatch.glerl.noaa.gov/glsea/glsea.html. To date, GLSEA2 serves as the best resource

243 to examine spatial and temporal variability of surface water temperature. Through the NOAA

244 CoastWatch program, GLSEA2 provides daily digital maps of the Great Lakes LST. The LST

245 data are stored as a 1024x1024 pixel map, available in PNG, ASC, and NetCDF formats, suitable

246 for visualization and further manipulation with readily available software. The data are derived

247 from NOAA/AVHRR (Advanced Very High Resolution Radiometer) satellite imagery and are

248 updated daily with information from the cloud-free portions of the previous day's satellite

249 imagery. To generate continuous evolution of the LST, a smoothing algorithm is applied to the

250 previous day's available map when no imagery are available, as described by Schwab et al

251 [1992]. In addition, direct in situ measurements of LST from the National Data Buoy Center

252 (NDBC, available from http://www.ndbc.noaa.gov/) are analyzed over the five Great Lakes to

11
253 assist in model evaluation. Moored water temperature data collected at the western basin of Lake

254 Superior for the year 2011, and maintained by University of Minnesota-Duluth [Titze and Austin,

255 2014], are also used to evaluate model simulations of subsurface water temperature for both

256 coupled and uncoupled 3-D circulation models. The mooring is located close to NDBC buoy

257 station 45006 at 47.335°N 89.793°W, at a local water depth of 183m.

258 Ice cover data over the Great Lakes are collected by the Canadian Ice Service and U.S.

259 National Ice Center (available from http://www.natice.noaa.gov/products/great_lakes.html) at a

260 temporal resolution of 3-7 days. The observations are constructed using various imagery sources,

261 with pixel resolutions down to 50 meters. The data also include necessary value-added

262 interpretation of these imagery sources to properly identify the extent of the ice edge contours.

263 Data used in the present study are linearly interpolated to daily time scales, and comparisons are

264 performed to evaluate model simulations of ice coverage and duration.

265 TC-3D-GLARM simulated precipitation and evaporation over the Great Lakes are

266 evaluated against lake-wide average estimates from the GLHCD dataset (only lake-average

267 monthly estimates are available), which assimilates nearshore, over-lake, and watershed-based

268 hydrometeorological measurements into model simulations of the major components of the water

269 budget for the Great Lakes basin [Hunter et al., 2015] . The GLHCD estimates over-lake

270 precipitation via spatial interpolation, using a modified version of Thiessen weighting and data

271 from stations that are located near shore or on islands and offshore lighthouses [Hunter et al.,

272 2015] . GLHCD monthly over-lake evaporation estimates for each lake are obtained from daily

273 simulations using NOAA-GLERL's 1-D Large Lake Thermodynamics Model (LLTM), forced by

274 aggregated nearshore and offshore hydrometeorological measurements [Hunter et al., 2015] .

275

12
276 3. Design of numerical simulations

277 RegCM4 simulations were conducted over North America for approximately three

278 decades from 1982-2014 (Table 2). We first ran RegCM4 as a standalone, uncoupled model for

279 the first two decades (1982-2001), initialized at 1 January 1982 using the first-year as spin-up.

280 During the period of standalone simulation (1982-2001), the SST of the ocean and LST of the

281 lakes were prescribed by weekly 0.5 NOAA daily optimum interpolation SST (OISST)

282 [Reynolds et al., 2007] . RegCM4 was then coupled with the hydrodynamic and ice models for

283 the remaining 13 years (2002-2014), using a 1-year spin-up (2002) for the hydrodynamic-ice

284 model. During the two-way coupled simulation period (2002-2014), the two model components

285 ran simultaneously, with coupled information exchange between RegCM4 and FVCOM at 3-

286 hour intervals. We note that tests with a higher coupling frequency showed no significant

287 improvement against more computational time. In the two-way coupled framework, the LST

288 fields and ice coverage are dynamically calculated by the 3-D hydrodynamic model and ice

289 model and provided to RegCM4 as over-lake surface boundary conditions. In turn, the surface

290 forcing fields required by the hydrodynamic and ice models are dynamically calculated and

291 provided by the coupled atmospheric model component. These fields include wind velocity,

292 precipitation, relative humidity, cloud cover, and incoming shortwave and longwave radiation.

293 Instantaneous estimates of latent and sensible heat flux and upward longwave radiation are made

294 within FVCOM at each time-step using simulated LSTs with the COARE 2.6 bulk algorithm

295 [Fairall et al., 1996] .

296 Although the MY25 sub-model is currently the most robust, commonly used turbulence

297 closure scheme in Great Lakes 3-D hydrodynamic models, the accurate simulation of vertical

298 mixing processes in water is still a challenging prospect in the hydrodynamic modeling

13
299 community. Over-mixing with the MY25 sub-model during the fall season has occasionally been

300 found to occur in the Great Lakes [Xue et al., 2015] . Hence, a nudging scheme [Ye et al. 2016]

301 is applied during November and December to constrain deep-water temperatures (> 100 m

302 depth) to 4 °C or cooler (the temperature of maximum density for fresh water), in case over-

303 mixing occurs during fall turnover.

304 There are several practical reasons why the above model configuration is used in order to

305 optimize the computational efforts: 1) For evaluating RegCM4’s performance in simulating

306 long-term climate trends and variability over North America, the full simulation time spanning

307 more than three decades (1982-2014) is clearly needed to provide an objective assessment; 2)

308 Due to computational loads, however, the TC-3D-GLARM simulation is performed only for

309 2003-2014 (i.e., 12 years). This is sufficient to appropriately evaluate the simulation of Great

310 Lakes hydrodynamic and ice thermodynamics in the new regional climate modeling system,

311 particularly given the large interannual variability during the 2012-2014 period; 3) As dimictic

312 lakes, the memory of initial conditions for the Great Lakes should generally be less than one

313 year; hence a 1-year spin-up (2002) of the hydrodynamic model is appropriate; and 4) high-

314 resolution ice atlas data (for model evaluation) are available from 2002 onward.

315 In order to assess the improvements of TC-3D-GLARM over other model configurations

316 (1-D coupling and uncoupled 3-D circulation models), we conducted three additional

317 experiments: one simulation with RegCM4 coupled with the 1-D lake model and two simulations

318 with the 3-D FVCOM model in a standalone, uncoupled configuration with RegCM4. Additional

319 details regarding the various simulation experiments are provided in Table 2.

320

321 4. Results

322 4.1 Simulated air temperature and precipitation

14
323 We begin by evaluating the modeled surface climatology of air temperature and

324 precipitation and assess the degree to which the model can capture regional climate trends and

325 variability, not only for the Great Lakes region but also across North America. The comparison

326 helps to provide confidence in the performance of the RCM, which is a prerequisite to ensure the

327 accurate representation of regional climate interactions with the Great Lakes hydrodynamic and

328 ice models.

329 As a measure of how well the model reproduces the observed interannual variability and

330 long-term trends in domain-averaged North American air temperature and precipitation, the

331 observed and simulated annual anomalies for both variables over the period 1982–2013 are

332 compared (Figure 3). The model accurately captures both the interannual variability and climatic

333 warming trend of surface air temperature over the past three decades. The simulated and

334 observed temperature anomalies are highly correlated, with a correlation coefficient of 0.95 and

335 a root-mean-square-error (RMSE) of 0.19 °C. The simulated long-term warming trend of 0.24 °C

336 per decade compares very well with the observed trend of 0.25 °C per decade, indicating that the

337 model captures the climatic trend as well as the interannual variability. In terms of annual

338 precipitation, the model similarly produces robust results in comparison with the CRU data, with

339 a high correlation coefficient of 0.84, a low RMSE of 24 mm (compared to an annual mean

340 precipitation of 748mm). Aside from a few exceptions (e.g., 2011), most of the wet years (e.g.

341 1983, 1990, 1993, 2004) and dry years (e.g. 1988, 2002, 2012) are well captured by the model.

342 Figure 4 presents a comparison of the simulated and observed spatial distribution of the

343 climatological annual mean surface air temperature and precipitation. The model reproduces the

344 observed spatial patterns in air temperature reasonably well, including large temperature

345 gradients over the Rocky Mountains in the western United States and a meander of warmer

15
346 temperature intrusions in the south and central United States (Figure 4a, b). The pattern

347 correlation between the modeled and observed climatology of surface air temperature is 0.83.

348 Differences between the modeled and observed temperature maps (Figure 4c) show close

349 agreement over the majority of North America, with errors of ~1° C throughout the Great Plains,

350 Midwest, Southwest, Southeast, and Northeast regions. Larger errors of up to 3 °C are present

351 over the Rocky Mountain region, which is likely due to the complex topography and model

352 resolution issues. Over the Great Lakes region, the modeled annual and seasonal (Figure S1) air

353 temperature climatology agrees very well with observations, with biases of less than 1 °C during

354 all seasons except winter (Figure S1). This exception is associated with a warm bias over the

355 northern boundary of the model domain during winter (DJF), which affects the water

356 temperature and ice cover simulation (section 4.2) of Lake Superior, particularly the northern

357 portion. A direct comparison of the spatial pattern of air temperature between the RegCM4

358 boundary input files generated from ERA-40 and the CRU observations suggests that the warm

359 bias in the northern boundary is most likely inherited from the driving ERA-40 reanalysis

360 product (not shown).

361 In terms of annual precipitation, both the spatial pattern and magnitude of the simulated

362 precipitation rate are consistent with observations (Figure 4d and 4e) with a pattern correlation of

363 0.88. This includes the general pattern of low precipitation to the west of the Missouri River,

364 higher precipitation to the east, and a narrow band of high precipitation along the Pacific

365 Northwest. The difference map (Figure 4f) reveals some wet biases in regions of topography and

366 a general dry bias in the southeastern U.S., but the Great Lakes region shows generally good

367 model-observation agreement, with most errors in seasonal precipitation being generally less

368 than 0.5-1.0 mm/day, aside from a few locations in lake-effect precipitation belts (Figure S2).

16
369

370 4.2 Great Lakes surface water temperature

371 As mentioned previously, our primary research goal in this study is to improve the

372 simulation of large lakes in regional climate modeling. Toward this end, we compare the model-

373 simulated and observed LST for each of the Great Lakes over the period 2003-2014 (Figures 5

374 and 6). At the climatological, mean-seasonal time scale, the modeled LST agrees well with

375 observations in each lake (Figure 5, left panels). For Lake Superior, the model shows a close

376 agreement with GLSEA2 with a RMSE of 0.88 °C. The most noticeable bias occurs during

377 summer (June-August). For Lake Michigan and Lake Huron, the difference in the mean seasonal

378 cycle of LST between the model and GLSEA2 are 0.65 oC and 0.79 oC, respectively, without a

379 particular bias in specific seasons. For Lake Ontario, there is a slight underestimate of summer

380 LST, and the RMSE is 0.93 °C. Relatively larger bias occurs in Lake Erie, where the model

381 shows an overestimate from June to December, with a RMSE of 1.43 °C. Reasons for such bias

382 are unclear and would require additional experiments to determine.

383 In addition to comparing the mean seasonal climatology of LST with observations for our

384 model assessment, we also present a model-data comparison of LST and ice cover anomalies

385 (deviations from the mean seasonal climatology) for 12 consecutive years (2003-2014) (Figures

386 6 and 9). This provides a more rigorous comparison, since it requires the model to reproduce

387 both the climatological mean seasonal cycle, as well as monthly to interannual variability. The

388 TC-3D-GLARM modeling system accurately captures the interannual variability of LST in all

389 five lakes, with much of the LST variability being significantly influenced by each lake’s depth

390 and geographic characteristics. Although the shallower lakes exhibit larger seasonal variability in

391 climatological LST (e.g. ~25 °C for Lake Erie compared to ~18°C for Lake Superior; Figure 5),

17
392 all five of the lakes exhibit similarly strong internnaual variability in LST, with a range of +/-

393 4°C (Figure 6). Such variability can become accentuated during extreme climatic events, such as

394 the cold winter of 2013-2014, which led to significantly reduced summer LST on Lake Superior

395 during 2014 (up to 4°C below its climatology value) due to significant delays in summer

396 stratification.

397 In addition to comparing the lake-wide mean LST with the remotely sensed GLSEA2

398 data, we also assessed TC-3D-GLARM’s performance against spatially distributed, in situ LST

399 measurements at all nine NDBC buoys on the five Great Lakes, each of which are located well

400 offshore, but in a variety of different water depths. For comparison, we also ran RegCM4 with

401 the default 1-D lake model. For brevity, we only show results for the last four years of the

402 simulation (Figure 7), which includes a very warm year (2012) and very cold year (2014). The 1-

403 D lake model distributed in the 2-D domain only capture the seasonal variability of LST at two

404 shallow sites – one in western Lake Erie (45005; 13m) and the other in southern Lake Huron

405 (45008; 54m). Otherwise, the 1-D model grossly fail to resolve the seasonality of LST at the

406 other seven stations with greater water depth. Similarly poor results from the default 1-D lake

407 model coupled in RegCM4 have been shown before in Figure 3 of Bennington et al., [2014] ,

408 who provided a very detailed assessment and further refinement of the 1-D lake model

409 simulation when coupled with RegCM4. In contrast, the new results from TC-3D-GLARM show

410 excellent agreement with each of the nine in situ buoy measurements on the Great Lakes,

411 regardless of water depth.

412 Due to the immense surface area and abruptly changing bathymetry, Great Lakes LSTs

413 vary significantly across each lake, particularly during the summertime. Figure 8 shows the

414 spatial pattern of the seasonal climatology of LST from TC-3D-GLARM and GLSEA2 data.

18
415 During the springtime (Figure 8a and 8b), the spatial variability of LST has just begun to develop

416 within each lake, and the temperature pattern reflects the impacts of variations in water depth and

417 latitude, resulting in earlier springtime warming in the southern lakes and shallower water. In

418 summer, strong thermal gradients continue to be evident across latitude and varying water depths

419 (e.g., coastal slope zones; Figure 8c and 8d). The spatial patterns of modeled summer LST are

420 particularly well captured for Lakes Huron-Michigan and Ontario, while the model slightly

421 overestimates LST in mid-lake portions of Lake Superior and southwestern portions of Lake

422 Erie.

423 During autumn, LST decreases to ~7 °C in Lake Superior and ~12 °C in Lake Erie due to

424 reductions in net radiation and increases in latent and sensible heat flux associated with stronger

425 winds and frequent passages of cold, dry air [Blanken et al., 2011; Xue et al., 2015] . The spatial

426 pattern of LST becomes fairly homogenous within each lake, except for Lake Michigan-Huron

427 due to its extensive orientation in the meridional direction (Figure 8e and 8f). Surface cooling

428 continues through winter (Figure 8g and 8h), accompanied by rapid ice formation in the

429 nearshore regions (discussed in the next section) and relatively lower ice coverage in mid-lake,

430 due both to larger heat content in deep water and to strong winds in the open water that can

431 retard ice formation [Assel, 1990; Wang et al., 2012] .

432

433 4.3 Ice simulation for the Great Lakes

434 In the Great Lakes, ice cover plays an important role in shaping the lakes’ energy and

435 water balance by affecting net radiation, evaporation (latent heat), and sensible heat flux.

436 Simulation of ice cover in the Great Lakes has long been a challenge, and there have been few

437 models (1-D or 3-D) that are able to address the issue satisfactorily [Dupont et al. 2012]. On one

438 hand, TC-3D-GLARM captures reasonably well the seasonal and interannual variability in ice

19
439 cover, as well as the differences among lakes (Figures 5, 9 and S4). On the other hand, peak areal

440 ice coverage in Lake Superior tends to be underestimated (Figure 5-f), likely due to the

441 aforementioned warm bias in winter air temperature over the Lake Superior region in RegCM4

442 (Figure S1). Both model and observations indicate that Lake Huron and Lake Michigan have

443 peak ice coverage of ~50% and ~25%, respectively, in the seasonal climatology (Figures 5-g, h),

444 with ice typically diminishing below ~10% by late April. Lake Ontario (Figure 5i) shows the

445 least ice cover among the five lakes – with a typical year having only 15-20% maximum areal

446 ice coverage – primarily because Lake Ontario is the second deepest of the five Great Lakes (by

447 mean depth), while also being in a warmer, southeastern location that is generally downwind of

448 the other lakes. As the shallowest lake, Lake Erie (Figure 5j) develops the highest ice coverage

449 (with mean peak coverage around 70-75%).

450 TC-3D-GLARM also captures the interannual variability of ice cover reasonably well,

451 showing large variability in all five lakes (Figure 9). For example, the anomalies in peak ice

452 coverage for Lake Superior are ~50-75% above normal during cold winters (2003, 2009, and

453 2014), while ice coverage may be ~20-30% below normal in warm years (2006, 2010, and

454 2012). Similar patterns are observed in other lakes except Lake Ontario, which shows lower

455 variability, similar to its lower mean values. Lake Erie often shows much larger intra-seasonal

456 variability, such as double-peak patterns during cold winters (2003 and 2014), abruptly changing

457 above- and below-normal ice cover in 2004, 2007, and 2008, and exceptional low-ice years

458 during 2006 and 2012. Interannual variability of ice cover on Lake Huron is also quite large,

459 with peak values that range from 60% above normal to 40% below normal, while the largest ice

460 cover anomalies on Lake Michigan are usually no more than 30% above normal, with 2014

20
461 (>90%) being the primary exception. Abnormally low ice cover during the warm years of 2006

462 and 2012 are also well captured by the model across all the lakes (Figure 9 and S4).

463 In addition to areal ice coverage, the duration of winter ice cover (from ice-on to ice-off,

464 which is defined with a threshold of 10% ice coverage within an observation pixel or a model

465 grid) is another key characteristic that is important to examine for the Great Lakes. The

466 climatology of observed ice duration from the National Ice Center for the period 2003-2014 is

467 compared with the TC-3D-GLARM simulation in Figure 10. The simulated spatial distribution

468 of ice duration agrees with observations in general, but some discrepancies exist in various

469 locations, particularly Lake Superior. Aside from the north shore of Minnesota, the model tends

470 to underestimate ice duration over the majority of Lake Superior, as would be expected from the

471 underestimated ice coverage (Figure 5) and aforementioned warm bias in winter air temperature.

472 The model does, however, capture the general pattern of longer ice duration in the shallow,

473 nearshore regions and shorter ice duration in the deeper, offshore locations – not only for Lake

474 Superior, but the other four lakes as well. This includes a number of shallow regions in Lake

475 Michigan-Huron with longer ice duration, such as Green Bay, Georgian Bay, the Straits of

476 Mackinac, and the southern shore of Lake Huron. The climatological pattern of ice duration in

477 Lake Ontario is also well captured (particularly the strong gradients in the northeastern portion of

478 the lake), as is the much higher ice duration in Lake Erie, and especially in the shallow, western

479 basin (~70-90 days).

480

481 4.4 Over-lake evaporation and precipitation

482 Over-lake evaporation and precipitation are important components of the energy and

483 water budgets of the Great Lakes, yet very few over-lake hydrometeorological measurements are

484 available [Blanken et al., 2011; Spence et al., 2013] . While this makes it extremely difficult to

21
485 derive measurement-based spatial patterns of over-lake precipitation and evaporation, it also

486 highlights the value of integrating observations and RCM modeling systems such as TC-3D-

487 GLARM for estimating Great Lakes surface energy and water fluxess [Hunter et al., 2015] .

488 Figure 11 presents a comparison of the simulated and observed 12-year monthly

489 climatology of mean over-lake precipitation and evaporation. Results from TC-3D-GLARM and

490 GLHCD show similar modeled and observed seasonality, with some variations among lakes that

491 reflect differences in atmospheric and hydrodynamic conditions. For example, the lowest

492 evaporation rates generally occur in early to late spring, when water temperatures are still cold,

493 but the air is becoming warmer and more humid [Lenters, 2004] . This minimum evaporation

494 occurs in April for Lake Erie, May for Lakes Michigan-Huron and Ontario, and June for Lake

495 Superior. The earlier minimum for Lake Erie occurs due to the fact that the lake warms and

496 stratifies much faster because of its shallow depth and southern location, which allows the lake-

497 air moisture gradient to evolve faster and in earlier seasons than the other lakes. Similarly, the

498 highest evaporation rates occur earlier for Lake Erie (October) and later for Lake Superior

499 (December-January), which reflects both the larger thermal inertia of Lake Superior and the

500 more extensive ice coverage on Lake Erie. Autumn Lake Erie evaporation rates in TC-3D-

501 GLARM are considerably lower than those from GLHCD LLTM, but it is suspected that this

502 reflects evaporation rate overestimate in the LLTM simulation, which are unrealistically higher

503 than even the highest rates on Lake Superior. Observed and simulated over-lake precipitation

504 generally display good agreement for all four lakes, including a slight tendency for higher

505 precipitation rates in late spring and early summer (e.g., Lakes Superior and Michigan-Huron).

506 Precipitation estimates from TC-3D-GLARM during the autumn seasons are generally lower

22
507 than the GLHCD, suggesting either a bias in the GLHCD spatial interpolation of land-based

508 stations or an underestimate of modeled precipitation.

509

510 5. Extreme climatic events: The “Big Heat” (2012) and “Big Chill” (2014)

511 During the winter of 2011-12, the United States experienced the fourth warmest winter

512 on record in more than a century (NOAA 2013). Located near the warming center, the Great

513 Lakes had the lowest ice coverage on record since the 1960’s (Figure 9 and S4), resulting in an

514 exceptionally early onset of stratification, a longer period of stratification, and a deeper

515 thermocline. Figure 12 (left panels) demonstrates the evolution of lake thermal structure during

516 2012 for each of the five Great Lakes, as simulated by the TC-3D-GLARM modeling system.

517 Separated by only one intervening year, the Great Lakes region then experienced an extremely

518 cold winter in 2013-14, caused by anomalous intrusions of Arctic air and high-amplitude wave

519 patterns in the jet stream (commonly referred to in the media as the polar vortex). As a result, the

520 Great Lakes experienced an extended period of ice coverage – even into the month of June for

521 Lake Superior – which was followed by a significantly delayed and shorter period of

522 stratification, along with a shallower thermocline (Figure 12, right panels).

523 These recent extreme climatic events have raised challenging questions for the regional

524 climate modeling community, as it is becoming apparent that strong interannual climatic

525 variability is inherent to the Great Lakes system (Van Cleave et al., 2014; Gronewold et al.,

526 2015). Thus, it is important to determine not only the extent to which Great Lakes climatic trends

527 (and impacts) can be predicted, but also to determine how well we can predict extreme

528 fluctuations.

529

23
530 6. Discussion

531 6.1 Necessity of resolving 3-D hydrodynamics

532 In many large-lake systems, including the Great Lakes, hydrodynamic processes that

533 control thermal structure are far more complicated than simple 1-D vertical mixing. The dynamic

534 regimes in the Great Lakes vary appreciably with changes in water depth and can be divided into

535 offshore water, the coastal boundary layer, and nearshore regions, similar to coastal oceans. Each

536 region is characterized by a different momentum balance that essentially determines the multi-

537 scaled flow structures in the Great Lakes. In the open water, the momentum balance is primarily

538 between the Coriolis force and barotropic (wind-driven) and baroclinic (thermal-driven) pressure

539 gradient forces, while in the coastal boundary layer, the bottom and lateral frictional forces play

540 an important role in the momentum balance.

541 As demonstrated in Figure 13, the mean horizontal flow structure in the near-surface

542 layer of the Great Lakes features multi-gyre circulations and strong coastal-jet currents with

543 strong spatial variability. Beletsky et al. [1999] and Rao and Schwab [2007] give excellent

544 summaries on these flow patterns and related physical processes. Furthermore, short-term flow

545 conditions and thermal structures during episodic events such as storms further complicate the

546 system, as transient surface velocities can reach >10-20 cm/s, with spatial temperature gradients

547 exceeding 0.01°C/m in the thermal bar region. In addition, the vertical flow structures that

548 directly impact and respond to mixing processes can also be extremely complex, such as double-

549 cell secondary circulation at thermal fronts [Chen et al. 2000] or two-layer baroclinic flow like

550 those that occur in the straits of Mackinac during the summer-stratified season [Anderson and

551 Schwab, 2016]. Also, surface inertia-gravity (Poincaré) waves, Kelvin waves, and topographic

552 waves can induce mixing across thermal gradients.

24
553 As mixing processes are closely associated with water velocity shear in both the

554 horizontal and vertical dimensions [Smagorinsky, 1963; Mellor and Yamada, 1982], such flow-

555 dependent mixing processes are unresolvable in a 1-D lake model, where mixing coefficients are

556 typically determined only as a function of wind. Due to the different spatial scales of

557 atmospheric and hydrodynamic processes and the complexities of hydrodynamic conditions, the

558 wind-based mixing coefficients are oversimplified, with much less accuracy in resolving spatial

559 variability in mixing processes. Therefore, although 1-D lake models may work well in small and

560 shallow lakes, where hydrodynamic processes are relatively homogeneous, large lake systems

561 almost always require substantial empirical calibration of eddy diffusivity with a factor that

562 could range from 102-105 at various water depth to compensate for not explicitly simulating 3-D

563 mixing processes that influence thermal transfer [Gu et al. 2015] or use a “virtual bottom cut-

564 off” often at 50 m to tackle this issue [Gula and Peltier, 2012; Subin et al., 2012 ].

565 Advective transport of heat is another equally important process for large lakes that is

566 unresolved in 1-D lake models. During the stratified season, significant wind events create

567 upwelling and downwelling of the thermocline due to Ekman transport. These circulation

568 patterns can create regions that transport significant heat. The wind-driven and density-driven

569 circulations also transport and redistribute heat within the lake. The resulting thermal gradients in

570 turn support and sustain circulation patterns and waves. These processes are completely

571 neglected in 1-D vertical thermal-balance models.

572 It is straightforward to see the importance of using 3-D, primitive equation, turbulent-

573 closure circulation models to represent the hydrodynamics (rather than 1-D lake models) across a

574 wide range of water depths in all five lakes, especially for simulating impacts on lake thermal

575 structure and subsequent lake-atmosphere interactions.

25
576

577 6.2 Inclusion of coupled lake-air dynamics in hydrodynamic simulations

578 Although 3-D circulation models have been fairly widely applied in “standalone” fashion

579 for a variety of Great Lakes and oceanographic applications, two major challenges are commonly

580 recognized by the hydrodynamic modeling community: 1) the uncertainty of surface forcing

581 representations and 2) the failure to effectively resolve air-sea feedback processes. A dynamic

582 representation of meteorological conditions has been suggested to reduce surface forcing

583 uncertainty in hydrodynamic modeling [Xue et al. 2015], which also suggests that resolving

584 water-air feedback processes should improve model performance.

585 To examine how much model improvement in the current study is due to the utilization

586 of a two-way coupled modeling approach, we conducted a “model assessment simulation” to

587 compare TC-3D-GLARM with an uncoupled 3-D hydrodynamic model simulation. We note that

588 we retained all the model configurations (i.e., model parameters and temporal and spatial

589 resolution) to be exactly the same as those in the original TC-3D-GLARM simulation.

590 To demonstrate the importance of two-way coupling for Great Lakes hydrodynamic

591 simulations, two scenarios are considered. In both scenarios we use the same regional climate

592 model and 3-D circulation model (FVCOM), but without the TC-3D-GLARM framework (i.e.,

593 uncoupled). The first scenario (scenario I) is intended for hydrodynamic forecasting/prediction

594 and does not include observed LSTs to drive the RCM. Rather, the RegCM model is first

595 coupled with default 1-D lake model to generate the surface forcing field that is then used to

596 drive the 3-D lake circulation model. The second scenario (scenario II) is a hydrodynamic

597 hindcast simulation, wherein the RegCM model was run with observed lake surface boundary

598 conditions from the GLSEA daily LST, and the output was then used to drive the 3-D FVCOM

26
599 model. In scenario II, the hydrodynamic simulation is expected to be better than scenario I, since

600 the output from the RegCM simulation has eliminated LST-induced error that may be present in

601 scenario I. Both cases were hot-started from coupled simulation results ending on 31 December

602 2009 and continue through 2011 (2-year simulation). We examine the models’ performance in

603 simulating the vertical structure of water temperature in Lake Superior at the western basin

604 offshore buoy location, and in comparison with observations and the TC-3D-GLARM

605 simulation.

606 Figure 14 presents a model-data comparison of the time evolution of the subsurface

607 thermal structure in the western basin of Lake Superior. The results demonstrate improvements

608 in the model simulation due to resolving lake-air feedbacks using a two-way coupled model.

609 Scenario I (panel b) shows the poorest simulation results, namely that the water stratifies much

610 earlier than suggested by the observations, and the mixing depth tends to deepen too rapidly from

611 summer into autumn, with warm water lingering in the subsurface until mid-December (i.e., later

612 than observed). In comparison, results from scenario II (panel c) show some improvement, with

613 more accurate onset of stratification and mixed-layer depth. The over-predicted warm water from

614 mid-November to December has also been corrected to a large degree. The most accurate

615 simulation (panel d) is provided by the TC-3D-GLARM model, showing a simulated thermal

616 structure that agrees well with observations. The model improves the accuracy of the onset of

617 stratification in mid-July. The warm bias in the upper 20m from September to October is

618 corrected, and the cooling of water after November is also more accurately simulated. The only

619 result that worsens slightly is the underestimate of near-surface temperature during the short

620 period from mid-October to early November. These comparisons not only validate the TC-3D-

621 GLARM simulation of subsurface water temperature, but also demonstrate the importance of

27
622 using a two-way coupled model approach for 3-D hydrodynamic simulation of the Great Lakes,

623 particularly as a tool for future prediction.

624 Another example that supports the integrated modeling approach is the reasonably

625 accurate simulation of ice cover and duration, despite the use of a relatively simple ice

626 thermodynamic model in the integrated modeling system. These results suggest that proper

627 simulation of complex ice structure on the Great Lakes is more strongly linked to the high

628 resolution, accuracy, and interactions between atmospheric conditions (e.g., surface heat fluxes

629 and winds) and hydrodynamic conditions (e.g., surface water temperature and currents) than to

630 the complexity of the ice model itself. In other words, without an accurate representation of both

631 the atmospheric and hydrodynamic components, it is implausible to deliver an accurate ice

632 simulation, regardless of the level of sophistication in the existing ice model.

633 The underlying mechanisms of how the representation of water-air interactions influence

634 both hydrodynamic and atmospheric model performances vary significantly depending on

635 characteristics of the regional climate [Wei et al., 2013; Turuncoglu et al., 2013; Song et al

636 2014; Xue et al. 2014; Xue et al. 2015]. Zhang et al. [1995], Terray et al. [2011] and Masson et

637 al. [2012] suggest short-term (e.g. intradaily, daily) variability of surface water temperature

638 could impact the regional climate by cascading short-term fluctuation into the larger scale

639 motion; Xue et al. [2014] and Xue and Eltahir [2015] demonstrate that simulations were

640 improved through direct local-scale air-sea feedbacks that primarily controlled vertical thermal

641 dynamic processes. Zhang et al. [1995] show that the improvement could be a combined

642 influence of local and large-scale processes. For large lakes, the local-scale feedbacks that

643 control over-lake heat flux components or heat redistribution are likely to play a dominant role

644 [Song et al. 2004 and Xue et al. 2015], rather than through feedbacks that modify large-scale

28
645 atmospheric circulation. This is also evident in our results through the comparison of the heat

646 budgets simulated by the TC-3D-GLARM and the RCM-coupled default 1-D lake models (next

647 section), since differences are primarily shown in variables that are directly affected by surface

648 water temperature.

649

650 6.3 Impact of 3-D hydrodynamic modeling on atmospheric variables

651 Comparing differences between the surface heat and water budgets simulated by the TC-

652 3D-GLARM system and the RCM-coupled default 1-D lake models (Figure S5-S9 and Table

653 S1), the results show that the most substantial differences are primarily for variables that are

654 directly affected by surface water temperature, including evaporation, sensible heat flux, and

655 upward longwave radiation. Though further improvements can certainly be made to the 1-D

656 models [Bennington et al. 2014], the unresolved hydrodynamics in the RCM-coupled 1-D model

657 are found to cause large biases in the LST simulation (e.g., Figure 6), which consequently causes

658 an excessive bias in annual evaporation amounts (e.g., an average difference of ~0.3 m across the

659 five Great Lakes, Table S1). The corresponding overestimate of annual latent heat flux is ~25

660 W/m2, compared to the estimates from TC-3D-GLARM. Sensible heat flux and upward

661 longwave radiation are also overestimated in the 1-D model by an average of 13 W/m2 and 5

662 W/m2, respectively. Noticeable but smaller differences of 0.17 mm/d are found in the

663 simulation of precipitation, which is likely due to the fact that both large-scale precipitation,

664 which is less impacted by change in LST and local convective precipitation, which is more

665 sensitive to local conditions contribute to the total precipitation. There are no significant changes

666 in other atmospheric variables such as solar radiation and cloud coverage when RegCM is

667 coupled with the 1-D lake model. We note that the above comparisons were made using seasonal

29
668 and annual time scales and basin-averaged values. It is anticipated that there should be more

669 noticeable impacts in the atmospheric variables at shorter time scales and finer spatial scales,

670 although such analyses are currently beyond the scope of this study.

671

672 7. Summary and Conclusions

673 In this study, we developed a regional lake-climate modeling system (TC-3D-GLARM)

674 comprised of a RCM coupled with a 3-D hydrodynamic lake and ice model, with a focus on the

675 long-standing issue of simulating large, deep lakes within a climate modeling setting. In our two-

676 way coupled modeling system, the water temperature, lake ice, and atmospheric surface

677 conditions are simulated simultaneously with different model components in TC-3D-GLARM,

678 and these variables are allowed to freely interact among each other. Results show that the

679 simulations of the Great Lakes thermal structure , surface fluxes, and ice are improved

680 significantly compared to previous studies that have utilized simpler modeling systems.

681 Characteristics of the atmosphere and the lake thermal structure, hydrodynamics, and ice are all

682 well captured by TC-3D-GLARM, indicating the unique capability of each modeling component

683 (as well as the integrated model) for accurately representing the various components of the

684 regional lake-climate system.

685 Despite these initial successes of the TC-3D-GLARM modeling system, there is still

686 plenty of room for future improvement. While the ice model adequately simulates the

687 thermodynamics of ice formation on the Great Lakes, horizontal transport and mechanical

688 redistribution of ice certainly play important roles in large lakes such as the Great Lakes. This is

689 even more crucial if the simulation is focused at short-time scales (e.g., hourly, daily, or weekly).

690 Future work will incorporate ice dynamics simulation into the coupled modeling system, once

30
691 the existing numerical code of ice dynamics is upgraded to make it more computationally

692 efficient. Model improvements can also be made in terms of the bulk parameterization of lake-air

693 heat, momentum, and moisture fluxes. Our ongoing research into different bulk

694 parameterizations suggests that TC-3D-GLARM currently produces reasonable estimates of

695 latent and sensible flux, but that these estimates may be particularly uncertain during late autumn

696 and winter (often when the fluxes are largest). Additional corrections, therefore, may provide

697 further improvements in the simulation of ice cover, water temperature, and hydrodynamics, not

698 only during cold months, but other seasons as well, particularly for the deeper lakes. In addition,

699 by converting runoff information from RegCM4 into river inflow in FVCOM, and building an

700 inter-connected five Great Lakes hydrodynamic model, TC-3D-GLARM will be able to provide

701 complete estimate of surface water cycle for the Great Lakes basin.

702

703 Acknowledgements

704 This is contribution 40 of the Great Lakes Research Center at Michigan Technological

705 University. Xue’s research was supported by the U.S. Environmental Protection Agency grant

706 GL00E01291-0 and the Michigan Tech Research Excellence Fund grant.

707

708

709

710

711

712

713

31
714 REFERENCES

715 Anderson, E. J. and D. J. Schwab (2013), Predicting the oscillating bi-directional exchange flow in the
716 Straits of Mackinac, J. Great Lakes Res., 39, 663-671.

717 Anderson, E.J., and D.J. Schwab (2016), Summertime Geostrophic Flow and Baroclinic Enhancement of
718 Freshwater Exchange in the Straits of Mackinac, JGR Oceans (in review)

719 Assel, A. A. (1990), An ice-cover climatology for lake erie and lake superior for the winter seasons 1897-
720 1898 to 1982-1983, Int. J. Climatol., 10, 731-748.

721 Austin, J. A. and S. M. Colman (2007), Lake Superior summer water temperatures are increasing more
722 rapidly than regional air temperatures: A positive ice-albedo feedback, Geophys. Res. Lett., 34, - L06604.

723 Bai, X., J. Wang, D. J. Schwab, Y. Yang, L. Luo, G. A. Leshkevich, and S. Liu (2013), Modeling 1993–
724 2008 climatology of seasonal general circulation and thermal structure in the Great Lakes using FVCOM,
725 Ocean Modelling, 65, 40-63.

726 Beletsky, D., D. Schwab, and M. McCormick (2006), Modeling the 1998-2003 summer circulation and
727 thermal structure in Lake Michigan, Journal of Geophysical Research: Oceans, 111, - C10010.

728 Bennington, V., M. Notaro, and K. D. Holman (2014), Improving Climate Sensitivity of Deep Lakes
729 within a Regional Climate Model and Its Impact on Simulated Climate, J. Climate, 27, 2886-2911.

730 Bitz, C. M. and W. H. Lipscomb (1999), An energy-conserving thermodynamic model of sea ice, Journal
731 of Geophysical Research: Oceans, 104, 15669-15677.

732 Blanken, P. D., C. Spence, N. Hedstrom, and J. D. Lenters (2011), Evaporation from Lake Superior: 1.
733 Physical controls and processes, J. Great Lakes Res., 37, 707-716.

734 Chen, C., R. C. Beardsley, and G. Cowles (2006), An unstructured grid, finite-volume coastal ocean
735 model (FVCOM) system, Oceanography, 19, 78-89.

736 Clites, A. H., J. Wang, K. B. Campbell, A. D. Gronewold, R. A. Assel, X. Bai, and G. A. Leshkevich
737 (2014a), Cold Water and High Ice Cover on Great Lakes in Spring 2014, Eos, Transactions American
738 Geophysical Union, 95, 305-306.

739 Clites, A. H., J. P. Smith, T. S. Hunter, and A. D. Gronewold (2014b), Visualizing relationships between
740 hydrology, climate, and water level fluctuations on Earth's largest system of lakes, J. Great Lakes Res.,
741 40, 807-811.

742 Dupont, F., P. Chittibabu, V. Fortin, Y. R. Rao, and Y. Lu (2012), Assessment of a NEMO-based
743 hydrodynamic modelling system for the Great Lakes, Water Qual. Res. J. Can., 47(3–4), 198–214,
744 doi:10.2166/Wqrjc.2012.014.

745 Fairall, C. W., E. F. Bradley, D. P. Rogers, J. B. Edson, and G. S. Young (1996), Bulk parameterization
746 of air-sea fluxes for Tropical Ocean-Global Atmosphere Coupled-Ocean Atmosphere Response
747 Experiment, Journal of Geophysical Research: Oceans, 101, 3747-3764.

32
748 Feser, F., B. Rockel, H. von Storch, J. Winterfeldt, and M. Zahn (2011), Regional Climate Models Add
749 Value to Global Model Data: A Review and Selected Examples, Bull. Amer. Meteor. Soc., 92, 1181-1192.

750 Fujisaki, A., J. Wang, H. Hu, D. J. Schwab, N. Hawley, and Y. R. Rao (2012), A modeling study of ice–
751 water processes for Lake Erie applying coupled ice-circulation models, J. Great Lakes Res., 38, 585-599.

752 Gao, G., C. Chen, J. Qi, and R. C. Beardsley (2011), An unstructured-grid, finite-volume sea ice model:
753 Development, validation, and application, Journal of Geophysical Research: Oceans, 116, - C00D04.

754 Giorgi, F. (2006), Regional climate modeling: Status and perspectives, J.Phys.IV France, 139, 101-118.

755 Giorgi, F. et al. (2012), RegCM4: model description and preliminary tests over multiple CORDEX
756 domains, Clim. Res., 52, 7-29.

757 Grell, G. A., J. Dudhia, and D. Stauffer (1994), A description of the fifth-generation Penn State/NCAR
758 Mesoscale Model (MM5). <br />, NCAR Technical Note, NCAR/TN-398+STR.

759 Grell, G. A. (1993), Prognostic Evaluation of Assumptions Used by Cumulus Parameterizations, Mon.
760 Wea. Rev., 121, 764-787.

761 Gronewold, A. D. et al. (2015), Impacts of extreme 2013-2014 winter conditions on Lake Michigan's fall
762 heat content, surface temperature, and evaporation, Geophys. Res. Lett., 42, - 2015GL063799.

763 Gronewold, A. D., A. H. Clites, J. P. Smith, and T. S. Hunter (2013), A dynamic graphical interface for
764 visualizing projected, measured, and reconstructed surface water elevations on the earth's largest lakes,
765 Environmental Modelling & Software, 49, 34-39.

766 Gu, H., J. Jin, Y. Wu, M. B. Ek, and Z. M. Subin (2015), Calibration and validation of lake surface
767 temperature simulations with the coupled WRF-lake model, Climatic Change, 375 129, 471-485, doi:
768 10.1007/s10584-013-0978-y.

769 Gula, J. and W. R. Peltier (2012), Dynamical Downscaling over the Great Lakes Basin of North America
770 Using the WRF Regional Climate Model: The Impact of the Great Lakes System on Regional Greenhouse
771 Warming, J. Climate, 25, 7723-7742.

772 Harris, I., P. D. Jones, T. J. Osborn, and D. H. Lister (2014), Updated high-resolution grids of monthly
773 climatic observations ? the CRU TS3.10 Dataset, Int. J. Climatol., 34, 623-642.

774 Holtslag, A. A. M., De Bruijn, E. I. F., and H. Pan (1990), A High Resolution Air Mass Transformation
775 Model for Short-Range Weather Forecasting, Mon. Wea. Rev., 118, 1561-1575.

776 Hostetler, S. W., G. T. Bates, and F. Giorgi (1993), Interactive coupling of a lake thermal model with a
777 regional climate model, Journal of Geophysical Research: Atmospheres, 98, 5045-5057.

778 Huang, A., Y. R. Rao, Y. Lu, and J. Zhao (2010), Hydrodynamic modeling of Lake Ontario: An
779 intercomparison of three models, Journal of Geophysical Research: Oceans, 115, - C12076.

780 Huang, A, Y. R. Rao, W. Zhang, (2012): On Recent Trends in Atmospheric and Limnological Variables
781 in Lake Ontario.J. Climate,25,5807-5816.

33
782 Hunter, T. S., A. H. Clites, K. B. Campbell, and A. D. Gronewold (2015), Development and application
783 of a North American Great Lakes hydrometeorological database — Part I: Precipitation, evaporation,
784 runoff, and air temperature, J. Great Lakes Res., 41, 65-77.

785 IPCC (2013), IPCC, 2013: climate change 2013: The physical science basis. contribution of working
786 group I to the fifth assessment report of the intergovernmental panel on climate change, in , edited by T.
787 F. Stocker, D. Qin, G. -K Plattner, M. Tignor, S. K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex, and P.
788 M. Midgley, pp. 1535, Cambridge University Press, Cambridge, United Kingdom and New York, NY,
789 USA.

790 Kiehl, J. T., J. J. Hack, G. B. Bonan, B. A. Boville, B. P. Briegleb, D. L. Williamson, and P. J. Rasch
791 (1996), Description of the NCAR Community Climate Model (CCM3)., NCAR Technical Note,
792 NCAR/TN-420+STR.

793 Lenters, J. D., J. B. Anderton, P. Blanken, C. Spence, and A. E. Suyker (2013), Assessing the Impacts of
794 Climate Variability and Change on Great Lakes Evaporation. In: 2011 Project Reports. D. Brown, D.
795 Bidwell, and L. Briley, eds.

796 Lenters, J. D. (2004), Trends in the Lake Superior Water Budget Since 1948: A Weakening Seasonal
797 Cycle, J. Great Lakes Res., 30, Supplement 1, 20-40.

798 Lipscomb, W. H. and E. C. Hunke (2004), Modeling Sea Ice Transport Using Incremental Remapping,
799 Mon. Wea. Rev., 132, 1341-1354.

800 Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki (2007), Ridging, strength, and stability in
801 high-resolution sea ice models, Journal of Geophysical Research: Oceans, 112, n/a-n/a.

802 Lofgren, B. M. (2014), Simulation of atmospheric and lake conditions in the Laurentian Great Lakes
803 region using the Coupled Hydrosphere-Atmosphere Research Model (CHARM)., NOAA Technical
804 Memorandum GLERL-165.

805 Lofgren, B. M. (2004), A model for simulation of the climate and hydrology of the Great Lakes basin,
806 Journal of Geophysical Research: Atmospheres, 109, - D18108.

807 Maykut, G. A. and N. Untersteiner (1971), Some results from a time-dependent thermodynamic model of
808 sea ice, Journal of Geophysical Research, 76, 1550-1575.

809 Mellor, G. L. and T. Yamada (1982), Development of a turbulence closure model for geophysical fluid
810 problems, Rev. Geophys., 20, 851-875.

811 NOAA National Centers for Environmental Information, State of the Climate: National Overview for
812 Annual 2012, published online January 2013.

813 Notaro, M., V. Bennington, and S. Vavrus (2015), Dynamically Downscaled Projections of Lake-Effect
814 Snow in the Great Lakes Basin, J. Climate, 28, 1661-1684.

815 Notaro, M., A. Zarrin, S. Vavrus, and V. Bennington (2013), Simulation of Heavy Lake-Effect
816 Snowstorms across the Great Lakes Basin by RegCM4: Synoptic Climatology and Variability*,+, Mon.
817 Wea. Rev., 141, 1990-2014.

34
818 O'Reilly, C. M. et al. (2015), Rapid and highly variable warming of lake surface waters around the globe,
819 Geophys. Res. Lett., - 2015GL066235.

820 Pal, J. S. et al. (2007), Regional Climate Modeling for the Developing World: The ICTP RegCM3 and
821 RegCNET, Bull. Amer. Meteor. Soc., 88, 1395-1409.

822 Pal, J. S., E. E. Small, and E. A. B. Eltahir (2000), Simulation of regional-scale water and energy
823 budgets: Representation of subgrid cloud and precipitation processes within RegCM, Journal of
824 Geophysical Research: Atmospheres, 105, 29579-29594.

825 Patterson, J. C. and P. F. Hamblin (1988), Thermal simulation of a lake with winter ice cover1, Limnol.
826 Oceanogr., 33, 323-338.

827 Y.R. Rao, D.J. Schwab(2007), Transport and mixing between the coastal and offshore waters in the Great
828 Lakes: a review J. Great Lakes Res., 33, pp. 202–218

829 Reynolds, R. W., T. M. Smith, C. Liu, D. B. Chelton, K. S. Casey, and M. G. Schlax (2007), Daily High-
830 Resolution-Blended Analyses for Sea Surface Temperature, J. Climate, 20, 5473-5496.

831 Rothrock, D. A. (1975), The energetics of the plastic deformation of pack ice by ridging, Journal of
832 Geophysical Research, 80, 4514-4519.

833 Schwab, D. J. and K. W. Bedford (1994), Initial Implementation of the Great Lakes Forecasting System:
834 A Real-Time System for Predicting Lake Circulation and Thermal Structure, Water Pollution Research
835 Journal of Canada, Water Poll. Res. J. Canada, 29, 203-220.

836 Schwab, D. J., G. A. Leshkevich, and G. C. Muhr (1992), Satellite Measurements of Surface Water
837 Temperature in the Great Lakes: Great Lakes Coastwatch, J. Great Lakes Res., 18, 247-258.

838 Shore, J. A. (2009), Modelling the circulation and exchange of Kingston Basin and Lake Ontario with
839 FVCOM, Ocean Modelling, 30, 106-114.

840 Smagorinsky, J. (1963), General circulation experiments with the primitive equations, Mon. Wea. Rev.,
841 91, 99-164.

842 Song, Y., F. H. M. Semazzi, L. Xie, and L. J. Ogallo (2004), A coupled regional climate model for the
843 Lake Victoria basin of East Africa, Int. J. Climatol., 24, 57-75.

844 Spence, C., P. D. Blanken, J. D. Lenters, and N. Hedstrom (2013), The Importance of Spring and
845 Autumn Atmospheric Conditions for the Evaporation Regime of Lake Superior, J. Hydrometeor, 14,
846 1647-1658.

847 Steiner, A., J. Pal, S. Rauscher, J. Bell, N. Diffenbaugh, A. Boone, L. Sloan, and F. Giorgi (2009), Land
848 surface coupling in regional climate simulations of the West African monsoon, Clim. Dyn., 33, 869-892.

849 Subin, Z. M., W. J. Riley, and D. Mironov (2012), An improved lake model for climate simulations:
850 Model structure, evaluation, and sensitivity analyses in CESM1, Journal of Advances in Modeling Earth
851 Systems, 4, - M02001.

35
852 Titze D. J., J. A. Austin (2014), Winter thermal structure of Lake Superior, Limnology and
853 Oceanography, 59, doi: 10.4319/lo.2014.59.4.1336.

854 Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony (1975), The thickness distribution of sea
855 ice, Journal of Geophysical Research, 80, 4501-4513.

856 Turuncoglu, U. U., G. Giuliani, N. Elguindi, and F. Giorgi (2013), Modelling the Caspian Sea and its
857 catchment area using a coupled regional atmosphere-ocean model (RegCM4-ROMS): model design and
858 preliminary results, Geosci. Model Dev., 6, 283-299.

859 Van Cleave, K., J. Lenters, J. Wang, and E. Verhamme (2014), A regime shift in Lake Superior ice
860 cover, evaporation, and water temperature following the warm El Niño winter of 1997-98, Limnology and
861 Oceanography, 59, 1889-1898.

862 Wang, G., M. Yu, J. S. Pal, R. Mei, G. B. Bonan, S. Levis, and P. E. Thornton (2015), On the
863 development of a coupled regional climate--vegetation model RCM--CLM--CN--DV and its validation in
864 Tropical Africa, Clim. Dyn., 46, 515-539.

865 Wang, J., X. Bai, H. Hu, A. Clites, M. Colton, and B. Lofgren (2012), Temporal and Spatial Variability
866 of Great Lakes Ice Cover, 1973-2010*, J. Climate, 25, 1318-1329.

867 Wang, J., H. Hu, D. Schwab, G. Leshkevich, D. Beletsky, N. Hawley, and A. Clites (2010),
868 Development of the Great Lakes Ice-circulation Model (GLIM): Application to Lake Erie in 2003–2004,
869 J. Great Lakes Res., 36, 425-436.

870 Wei, J., P. Malanotte-Rizzoli, E. B. Eltahir, P. Xue, and D. Xu (2013), Coupling of a regional
871 atmospheric model (RegCM3) and a regional oceanic model (FVCOM) over the maritime continent,
872 Clim. Dyn., 43, 1575-1594.

873 Williamson, J. (1854), The Inland Seas of North America ; and, the Natural and Industrial Productions of
874 Canada with the Real Foundations for its Future Prosperity by James Williamson., H. Ramsay, Kingston
875 [Ont.] :.

876 Xue, P., C. Chen, and R. C. Beardsley (2012), Observing system simulation experiments of dissolved
877 oxygen monitoring in Massachusetts Bay, Journal of Geophysical Research: Oceans, 117, - C05014.

878 Xue, P., C. Chen, R. C. Beardsley, and R. Limeburner (2011), Observing system simulation experiments
879 with ensemble Kalman filters in Nantucket Sound, Massachusetts, Journal of Geophysical Research:
880 Oceans, 116, - C01011.

881 Xue, P. and E. A. B. Eltahir (2015), Estimation of the Heat and Water Budgets of the Persian (Arabian)
882 Gulf Using a Regional Climate Model, J. Climate, 28, 5041-5062.

883 Xue, P., E. A. B. Eltahir, P. Malanotte-Rizzoli, and J. Wei (2014), Local feedback mechanisms of the
884 shallow water region around the Maritime Continent, Journal of Geophysical Research: Oceans, 119,
885 6933-6951.

36
886 Xue, P., D. J. Schwab, and S. Hu (2015), An investigation of the thermal response to meteorological
887 forcing in a hydrodynamic model of Lake Superior, Journal of Geophysical Research: Oceans, 120,
888 5233-5253.

889 Ye. X. and P. Xue. (2016), Impacts of vertical mixing on the lake thermal structure and ice formation in a
890 hydrodynamic model. Journal of coastal research (submitted)

891 Zeng, X., M. Zhao, and R. E. Dickinson (1998), Intercomparison of Bulk Aerodynamic Algorithms for
892 the Computation of Sea Surface Fluxes Using TOGA COARE and TAO Data, J. Climate, 11, 2628-2644.

893 Zhong, Y., Notaro, M., Vavrus, S. J. and Foster, M. J. (2016), Recent accelerated warming of the Laurentian Great
894 Lakes: Physical drivers. Limnol. Oceanogr.. doi:10.1002/lno.10331

895

896

897

898

899

900

901

902

903

904

905

906

907

908

909

910

911

912

37
913 Table 1. Morphometric information for the Great Lakes.

Lake Lake Lake Lake Lake

Superior Michigan Huron Ontario Erie
Average Depth (m) 149 85 59 86 19

Maximum Depth (m) 406 281 229 244 64

Volume (km3) 12,232 4,918 3,538 1,639 483

Water Surface Area (km2) 82,097 57,753 59,565 19,009 25,655

914

915

916

917

918

919

920

921

922

923

924

925

926

927

928

929

930

38
931 Table 2. Summary of numerical model simulations.

Simulation Meteorological
Model LST boundary
Period (spin- forcing of lake Purpose
Configuration condition
up excluded) models

RegCM4 20 years Prescribed by Part of RCM

N/A
standalone (1983-2002) daily OISST validation

Simulated by the Simulated by the Part of RCM

12 years coupled 3-D coupled RCM with validation;
TC-3D-GLARM
(2003-2014) circulation 3-D TC-3D-GLARM
model hydrodynamics validation

Simulated by the Intercomparison

Simulated by the
RegCM coupled 4 years coupled RCM with of RCM coupled
coupled 1-D lake
with 1-D lake (2011-2014) 1-D thermal with 1-D vs 3-D
model
dynamics lake model
Intercomparison
3-D lake Simulated by of coupled vs
circulation model 2 years RegCM4 coupled uncoupled 3-D
N/A
(FVCOM) (2010-2011) with 1-D lake circulation model
standalone model (forecast
scenario)
Intercomparison
3-D lake Simulated by of coupled vs
circulation model 2 years RegCM4 with LST uncoupled 3-D
N/A
(FVCOM) (2010-2011) prescribed by circulation model
standalone GLSEA (hindcast
scenario)
932

933

934

935

936 Figure Captions

39
937 Figure 1: Upper panel: North American RCM model domain with topography. Lower panel:

938 Bathymetry of the North American Great Lakes with NDBC buoy station locations

939 (black squares).

940 Figure 2: Unstructured triangular mesh used in the Great Lakes 3-D hydrodynamic model.

941 Figure 3: Annual temperature anomalies (upper panel) and precipitation anomalies (lower panel)

942 over North America for the period 1983–2013 from the model simulation and CRU

943 observations.

944 Figure 4: Model-data comparison showing the climatology of surface air temperature and

945 precipitation over North America (1983-2013). Left panel: CRU observations; Middle

946 panel: model simulation; Right panel: Difference (model minus observations).

947 Figure 5: Climatology (2003-2014) of daily LST and ice coverage from the TC-3D-GLARM

948 simulation (red dashed lines) and observations (blue lines; LST from GLSEA2 and ice

949 cover from the National Ice Center).

950 Figure 6: Time series of lake-average LST anomalies (monthly value minus long-term monthly

951 mean; 2003-2014), as simulated by TC-3D-GLARM (red) and compared with

952 GLSEA2 observations (blue).

953 Figure 7: Time series of surface water temperature simulated by TC-3D-GLARM (red) and

954 RegCM coupled with 1-D lake models (black), in comparison with in situ

955 measurements at 9 NDBC buoy stations (blue; see Figure 1 for buoy locations) for

956 2011-2014.

957 Figure 8: Model-data comparison of the LST seasonal climatology (2003-2014) during spring

958 months (AMJ), summer months (JAS), fall months (OND), and winter months (JFM).

959 Left panel: GLSEA2 observations; Right panel: TC-3D-GLARM simulation.

40
960 Figure 9: Time series of lake-average ice cover anomalies (daily value minus long-term daily

961 mean; 2003-2014), as simulated by TC-3D-GLARM (red) and compared with

962 observations (blue) from the National Ice Center.

963 Figure 10: Climatology of observed ice duration on the Great Lakes (in days) from the National

964 Ice Center (top) and the TC-3D-GLARM simulation (middle) for the period 2003-

965 2014. Bottom panel shows the difference (model minus observed).

966 Figure 11: Seasonal climatology of monthly mean over-lake precipitation (red) and evaporation

967 rate (blue) from GLHCD observations and TC-3D-GLARM simulation results (2003-

968 2014).

969 Figure 12: Simulated Great Lakes water temperature profiles from the TC-3D-GLARM model

970 at various buoy locations (Figure 1) on Lake Superior (45001), Lake Michigan

971 (45007), Lake Huron (45008), Lake Ontario (45012), and Lake Erie (45007) during

972 the “Big Heat” of 2012 (left) and “Big Chill” of 2014 (right).

973 Figure 13: Summer-mean (JAS)-mean circulation patterns on the Great Lakes (upper 20 m), as

974 simulated by the TC-3D-GLARM model for the period 2003-2014.

975 Figure 14: Model-data comparison of time evolution of the lake subsurface thermal structure in

976 the western basin of Lake Superior for year 2011 (see section 2 for mooring data

977 description).

978

979

41
Figure 1 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 1.pdf

48 o 4900
N 3000
1800
1100
670
40 o 400
N 250

(meter)
150
90
32 o 55
N 33
20
12
7
24 o 4
N 3
2
1

oW
135 o 60
W
o
120 oW 75 W
o o
105 W 90 W

Fig. 1 Upper panel: North American RegCM model domain with topography. Lower panel: Bathymetry of the
North American Great Lakes with NDBC buoy station locations (black squares).

1
Figure 2 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 2.pdf

49

48

47

46
Latitude (degree)

45

44

43

42

41
-92 -90 -88 -86 -84 -82 -80 -78 -76
Longitude (degree)

Fig. 2 Unstructured triangular mesh used in the Great Lakes 3-D hydrodynamic model.

2
Figure 3 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 3.pdf

Surface Air Temperature Anomaly(oC) 2

corr= 0.95
1 rmse= 0.19 oC

−1
RegCM
CRU
−2
1983 1988 1993 1998 2003 2008 2013

200
Precipitation Anomaly (mm/year)

150 corr= 0.84

100 rmse= 24.19 mm/year
50

−50

−100
RegCM
−150
CRU
−200
1983 1988 1993 1998 2003 2008 2013
Time(year)

Fig. 3 Annual temperature anomalies (upper panel) and precipitation anomalies (lower panel) over North
America for the period 1983–2013 from the model simulation and CRU observations.

3
Figure 4 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 4.pdf

CRU 1983−2013 RegCM 1983−2013

55

50

45
LATITUDE

40

35

30

25
(a) (b) (c)
20
−130 −120 −110 −100 −90 −80 −70 −130 −120 −110 −100 −90 −80 −70
LONGITUDE LONGITUDE LONGITUDE

T(oC)
−8 −4 0 4 8 12 16 20 24 28

(d) (e) (f)

Fig. 4 Model-data comparison showing the climatology of surface air temperature and precipitation over North
America (1983-2013). Left panel: CRU observations; Middle panel: model simulation; Right panel: Difference
(model minus observations).

4
Figure 5 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 5.pdf

Climatology LST Climatology Ice Cover

80
25 GLSEA2 Lake Superior
o 70
Observation Lake Superior
TC−3D−GLARM rmse= 0.88 C TC−3D−GLARM rmse= 2.29 %
60
20
(f)

Ice Cover(%)
(a) 50
T(oC)

15
40

10 30

20
5
10

0 0
80
Lake Michigan Lake Michigan
25
(b) o
rmse= 0.65 C 70 (g) rmse= 1.75 %
60
20

Ice Cover(%)
50
T(oC)

15
40

10 30

20
5
10

0 0
80
Lake Huron Lake Huron
25
(c) o
rmse= 0.79 C 70 (h) rmse= 2.49 %
60
20
Ice Cover(%)

50
T(oC)

15
40

10 30

20
5
10

0 0
80
Lake Ontario Lake Ontario
(i)
25
(d) o
rmse= 0.93 C 70
rmse= 2.66 %
60
20
Ice Cover(%)

50
T(oC)

15
40

10 30

20
5
10

0 0
80
Lake Erie o Lake Erie
25
(e) rmse= 1.43 C 70 (j) rmse= 4.47 %
60
20
Ice Cover(%)

50
T(oC)

15
40

10 30

20
5
10

0 0
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

Fig. 5 Climatology (2003-2014) of daily LST and ice coverage from the TC-3D-GLARM simulation (red
dashed lines) and observations (blue lines; LST from GLSEA2 and ice cover from the National Ice
Center).

5
Figure 6 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 6.pdf

Temperature Anomaly (oC) 6

GLSEA2 Lake Superior

4 TC−3D−GLARM

−2

−4
4
Temperature Anomaly (oC)

3 Lake Michigan
2

−1

−2

−3

−4
4
Temperature Anomaly (oC)

3 Lake Huron
2

−1

−2

−3

−4
4
Temperature Anomaly (oC)

3 Lake Ontario
2

−1

−2

−3

−4
4
Temperature Anomaly (oC)

3 Lake Erie
2

−1

−2

−3

−4
2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014
Time(year)
Fig. 6 Time series of lake-average LST anomalies (monthly value minus long-term monthly mean; 2003-2014),
as simulated by TC-3D-GLARM (red) and compared with GLSEA2 observations (blue).
6
Figure 7 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 7.pdf

Buoy
1−D−Coupled
TC-3D-GLARM
30
Lake Superior 45001 (255m) Lake Superior 45006 (195m) Lake Superior 45004 (238m)
25

20
T(oC)

15

10

0
30
Lake Michigan 45002 (181m) Lake Michigan 45007 (159m) Lake Ontario 45012 (145m)
25

20
T(oC)

15

10

0
30
Lake Huron 45003 (135m) Lake Huron 45008 (54m) Lake Erie 45005 (13m)
25

20
T(oC)

15

10

0
2011 2012 2013 2014 2011 2012 2013 2014 2011 2012 2013 2014
Time(year) Time(year) Time(year)

Fig. 7 Time series of surface water temperature simulated by TC-3D-GLARM (red) and RegCM coupled with 1-
D lake models (black), in comparison with in-situ measurements at 9 NDBC buoy stations (blue; see Figure 1
for buoy locations) for 2011-2014.

7
Figure 8 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 8.pdf

GLSEA2 LST TC-3D-GLARM

(a) (b)

Latitude (degree)
48

46

44

42
AMJ

(c) (d)
Latitude (degree)

Latitude (degree)
48

46

44

42
JAS

(e) (f)
Latitude (degree)

48

46

44

42
OND

(g) (h)
Latitude (degree)

Latitude (degree)

48

46

44

42
JFM
−90 −85 −80 −90 −85 −80
Longitude (degree) Longitude (degree)

T(oC)
0 5 10 15 20 25 30

Fig. 8 Model-data comparison of the LST seasonal climatology (2003-2014) during spring months (AMJ),
summer months (JAS), fall months (OND), and winter months (JFM). Left panel: GLSEA2 observations; Right
panel: TC-3D-GLARM simulation.

8
Figure 9 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 9.pdf

80
Observation
Ice Cover Anomaly(%)

60
TC−3D−GLARM Lake Superior
40

20

−20

−40

−60

80
Ice Cover Anomaly(%)

60 Lake Michigan
40

20

−20

−40

−60

80
Ice Cover Anomaly(%)

60 Lake Huron
40

20

−20

−40

−60

80
Ice Cover Anomaly(%)

60 Lake Ontario
40

20

−20

−40

−60

80
Ice Cover Anomaly(%)

60 Lake Erie
40

20

−20

−40

−60

2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014
Time(year)
Fig. 9 Time series of lake-average ice cover anomalies (daily value minus long-term daily mean; 2003-2014), as
simulated by TC-3D-GLARM (red) and compared with observations (blue) from the National Ice Center.
9
Figure 10 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 10.pdf

Climatology of Ice Duration

Fig. 10 Climatology of observed ice duration on the Great Lakes (in days) from the National Ice Center (top) and the
TC-3D-GLARM simulation (middle) for the period 2003-2014. Bottom panel shows the difference (model minus
observed).

10
Figure 11 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 11.pdf

3
mm/d

1
TC−3D−GLARM P
GLHCD P
0 TC−3D−GLARM E
GLHCD E
Lake Superior
−1
5

3
mm/d

Lake Michigan−Huron
0
5

3
mm/d

Lake Ontario
0
8

5
mm/d

0
Lake Erie
−1
Jan Feb March April May June July August Sep Oct Nov Dec

Fig. 11 Seasonal climatology of monthly mean over-lake precipitation (red) and evaporation rate
(blue) from GLHCD observations and TC-3D-GLARM simulation results (2003-2014).

11
Figure 12 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 12.pdf

2012 Temperature Profile 2014 Temperature Profile

Depth(m)−Lake Superior

−50 −50

−100 −100

−150 −150

−200 −200
Depth(m)−Lake Michigan

−20 −20

−40 −40

−60 −60

−80 −80

−100 −100

−120 −120

−140 −140

−5 −5
Depth(m)−Lake Huron

−10 −10
−15 −15
−20 −20
−25 −25
−30 −30
−35 −35
−40 −40
−45 −45
Depth(m)−Lake Ontario

−20 −20

−40 −40

−60 −60

−80 −80

−100 −100

−120 −120

−2 −2
Depth(m)−Lake Erie

−4 −4

−6 −6

−8 −8

−10 −10

−12 −12

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

o
T( C)
0 5 10 15 20 25

Fig. 12 Simulated Great Lakes water temperature profiles from the TC-3D-GLARM model at various buoy locations
(Figure 1) on Lake Superior (45001), Lake Michigan (45007), Lake Huron (45008), Lake Ontario (45012), and Lake
Erie (45007) during the “Big Heat” of 2012 (left) and “Big Chill” of 2014 (right).

12
Figure 13 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 13.pdf

49

Summer Mean Circulation

0.1 m/s
48

47

46
Latitude (degree)

45

44

43

42

41
-92 -90 -88 -86 -84 -82 -80 -78 -76
Longitude (degree)

Fig. 13 Summer (JAS)-mean circulation patterns on the Great Lakes (upper 20 m), as simulated by the TC-3D-
GLARM model for the period 2003-2014.

13
Figure 14 Click here to download Rendered Figure Figures_TC-
GLARM_revision_xy 14.pdf

Observation Scenario I
0 0
4
−10 −10
Depth(m)

−20 −20

−30 −30

−40 −40

−50 −50

Scenario II TC-3D-GLARM
0 0

−10 −10
Depth(m)

−20 −20

−30 −30

−40 −40

−50 −50
Jun/15 Jul/15 Aug/15 Sep/15 Oct/15 Nov/15 Dec/15 Jun/15 Jul/15 Aug/15 Sep/15 Oct/15 Nov/15 Dec/15
Time Time

T(oC)
5 8.75 12.5 16.25 20

Fig. 14 Model-data comparison of time evolution of the lake subsurface thermal structure in the western basin of
Lake Superior for year 2011 (see section 2 for mooring data description).

14