Developing the Coupled CWRF‐FVCOM Modeling System to Understand and Predict Atmosphere‐Watershed Interactions Over the Great Lakes Region

Abstract

Coupling 3‐D hydrodynamics with climate models is necessary but difficult for resolving multiscale interactions and has been rarely implemented in predicting Great Lakes' water level fluctuations because of issues in treating net basin supply (NBS) components and connecting channel flows. This study developed an interactive lake‐atmosphere‐hydrology modeling system by coupling the regional Climate‐Weather Research and Forecasting model (CWRF) with the 3‐D unstructured‐grid Finite Volume Coastal Ocean Model (FVCOM) in the Great Lakes region. The sensitivity of the coupled system, relative to the CWRF baseline using the 1‐D Lake, Ice, Snow and Sediment Simulator (LISSS), was evaluated in representing lake‐climate conditions during 1999–2015 against observations. As coupled with CWRF, FVCOM outperformed LISSS in simulating water surface temperature, ice cover, and vertical thermal structure at seasonal to interannual scales for all the five lakes and realistically reproduced the regional circulation patterns. In warm seasons, the improved lake conditions significantly corrected LISSS overestimates of surface air temperature together with larger‐scale circulation changes. Consequently, precipitation was generally reduced over each lake basin, mainly because of decreased surface moisture and heat fluxes along with enhanced atmospheric stability. Through the dynamic coupling, FVCOM predicts the water level fluctuations in direct response to the CWRF NBS components and connecting channel flows based on a stage‐fall‐discharge formulation. This coupled CWRF‐FVCOM reasonably captured the NBS variations and predicted the water level fluctuations for Lakes Superior and Michigan‐Huron. It represents a major advance in interacting regional climate and watershed processes to dynamically predict Great Lakes' water level seasonal‐interannual variations.

1 Introduction

The Great Lakes, consisting of Lakes Superior, Michigan, Huron, Erie, and Ontario with a combined surface area of approximately 246,000 km², are the largest fresh surface water source in the world. Their tremendous horizontal scales, characterized by distinct radiative and thermal properties compared to surrounding lands, considerably affect the local weather and regional climate through surface-atmospheric interactions (Bennington et al., 2014; Notaro et al., 2013; Scott & Huff, 1996). Due to lake's huge heat capacities, air temperature's spring warming and autumn cooling in the Great Lakes region are delayed with stifled variations both at diurnal and seasonal scales (Bates et al., 1993; Notaro et al., 2013). Significant air-water temperature differences in late fall and winter increase the atmospheric instability and thus induce the lake-effect precipitation on the lee side (Niziol et al., 1995; Scott & Huff, 1996). The lake-effect precipitation would be weakened in high ice-cover winter due to the reduced turbulent heat and moisture fluxes (Gerbush et al., 2008; Vavrus et al., 2013). Moreover, the socioeconomic and environmental effects of changing lake water levels over the coastal regions have been well recognized, ranging from coastal erosion and infrastructure damage by high water levels, to the constraint on shipping, port activities, and tourism during low-water periods (Dawson & Scott, 2010; Millerd, 2011; Wall, 1998).

Climate variability and change in turn significant influence lakes' thermal structure (Austin & Colman, 2007, 2008; King et al., 1997; Van Cleave et al., 2014) and water levels (Hanrahan et al., 2010). Over the past decades, under the general trend of global warming, an accelerated warming of water temperature was observed (Zhong et al., 2016). The resulting decline of ice cover could strengthen the heating process through the positive ice-albedo feedback (Austin & Colman, 2007). Meanwhile, the Great Lakes water level has been declining since the late 1980s, along with enhanced variability due to increased precipitation extremes (Gronewold & Stow, 2014) and is expected to continue declining in the future (Angel & Kunkel, 2010; Lofgren & Rouhana, 2016). Consequently, an integrated regional modeling system is of great importance and urgency to improve the understanding of lake-land-atmosphere interactions and water balance over the Great Lakes region, predict lake thermodynamic variations and water level fluctuations under a changing climate, and assess the climate change impacts on environment, economy, human health, and social activities at local and regional scale (Sharma et al., 2018).

In the past two decades, there have been a growing number of attempts to develop regional coupled atmosphere-ocean models for both weather and climate research (Li et al., 2014; Peng et al., 2012; Seo et al., 2007; Turuncoglu et al., 2013; Wei et al., 2014). Recently, the popular Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) has been coupled with various ocean models and implemented in several regions of the globe oceans (Masson et al., 2012; Samala et al., 2013; Warner et al., 2010; Zhao & Chan, 2017). In principle, these regional coupled models can be used to represent lake-atmosphere interactions over large, deep lakes (Anyah et al., 2006; Long et al., 2007; Song et al., 2004; X. Sun et al., 2014). However, there are important limitations for doing so. First, almost all these oceanic models were designed on regular grids and coupled with atmospheric models tightly on the same grids or through flux couplers onto coarser grids. The use of a uniform grid makes it extremely difficult to effectively resolve the wide range of motions that are essential to water thermodynamic structures and lake-atmospheric interactions. It is more challenging for such models to capture the coastal processes where lakes are interacting with not only the overlaying atmosphere but also the neighboring watersheds with complex and critical lateral boundary structures (such as runoff and channel flows). Second, most of the atmospheric models were outdated or lacking key watershed processes (Liang et al., 2019). Even the modern WRF, which was originally designed for short-range numeric weather prediction, requires substantial improvements in physics representations for long-term climate applications (Liang et al., 2012). Third, none of these coupled models was designed with a rigorous consideration of lateral interactions with terrestrial hydrology, which are essential for water level prediction.

To date, a number of one-dimensional (1-D) lake models, including two-layer Freshwater Lake model (Mironov et al., 2010) based on self-similarity theory and multilevel thermal models like Hostetler model (Hostetler et al., 1993) and Lake, Ice, Snow and Sediment Simulator (LISSS, Subin et al., 2012) based on wind-driven eddy diffusion, have been developed and coupled with regional climate models to investigate the atmosphere-lake-ice interactions in the Great Lakes region (Gula & Peltier, 2012; Notaro et al., 2013; Vavrus et al., 2013; Wright et al., 2013). In spite of a series of modifications and improvements (Bennington et al., 2014; Gu et al., 2015; Xiao et al., 2016), 1-D lake models still show obvious shortcomings in reproducing the lake thermal structure especially for large, deep lakes due to the absence of lake horizontal processes and hydrodynamics (Bennington et al., 2014; Mallard et al., 2014; Song et al., 2004). Meanwhile, a set of 3-D hydrodynamic models, including Princeton Ocean Model (Beletsky et al., 2006, 2013; Fujisaki et al., 2013; Huang et al., 2010), Nucleus of European Modeling of the Ocean (Dupont et al., 2012), and Finite Volume Coastal Ocean Model (FVCOM, Bai et al., 2013; C. Chen et al., 2003; Mao & Xia, 2017; Nguyen et al., 2017; Niu et al., 2015; Xue et al., 2015), have been performed in one-way off-line simulations under prescribed meteorological forcings for one or more Great Lakes with reasonable success. Xue et al. (2017) were the first to develop the two-way coupling of the Regional Climate Modeling system Version 4 (RegCM4) with FVCOM over the entire Great Lakes, which improved simulating lake thermal conditions over a 1-D lake model. However, the altered lake effects on regional climate, compared to the 1-D lake coupling, were only evaluated with a 4-year simulation, insufficient for climate variability. More importantly, the coupled system was unable to predict the water level fluctuations of the Great Lakes because it did not account for net basin supply (NBS) components (including overlake precipitation and evaporation, and lateral land surface runoff) and connecting channel flows. That is, the land-lake hydrology interaction was missing. Durnford et al. (2018) pioneered in developing the operational water cycle prediction system linking atmospheric, surface, river, and lake components over the Great Lakes and St. Lawrence River watershed. However, this system operates at a short to medium range and has not been evaluated for its suitable application on climate scales.

The Climate extension of WRF (CWRF, Liang et al., 2012) not only significantly improves regional climate simulations over RegCM4.6 (Liang et al., 2019) but also has built-in a unique Conjunctive Surface-Subsurface Process (CSSP) model that provide lateral flows from land to lakes to facilitate the watershed runoff contribution to NBS. The primary goal of this paper is to develop an interactive lake-atmosphere-hydrology modeling system in the Great Lakes region by coupling CWRF (Liang et al., 2012) with FVCOM (C. Chen et al., 2003) through the CPL7 coupler (Craig et al., 2012) and use it to understand how such dynamic coupling affects regional climate and lake water predictions. As a unique contribution to the literature, FVCOM is also constructed to dynamically predict the water level changes based on the NBS components provided by CWRF and connecting channel flows calculated using a stage-fall-discharge formulation. We then evaluate the ability of the coupled system in representing water surface temperature, ice cover, vertical thermal structure, and circulation patterns relative to the CWRF baseline with built-in LISSS (Subin et al., 2012) against available observations, as well as the regional climate responses during 1998–2015. Finally, we examine the system ability to capture the Great Lakes' water level variations and discuss future refinements.

2 Model System Components and Coupling 2.1 Atmosphere and Surface Interactive Component (CWRF)

CWRF has been continuously developed as a Climate extension of WRF (Skamarock et al., 2008) with numerous advances in the representation of essential physical processes, including interactions between land-atmosphere-ocean, convection-microphysics, cloud-aerosol-radiation (CAR), and system consistency throughout all process modules (Liang et al., 2012). It has been demonstrated with outstanding performance in capturing observed surface radiation, terrestrial hydrology, precipitation and extremes, and regional climate predictions (L. Chen et al., 2016; Liang et al., 2012, 2019; Liu et al., 2015; C. Sun & Liang, 2020; Yuan & Liang, 2011b). For the atmosphere, CWRF has a built-in CAR ensemble model that incorporates a wide variety of alternate parameterizations for cloud and aerosol properties, radiation transfers, and their interactions (Liang & Zhang, 2013; Zhang et al., 2013). It also incorporates an ensemble cumulus parameterization (ECP, Qiao & Liang, 2015, 2016, 2017), rooted in the G3 scheme (Grell & Dévényi, 2002) with updates in closure assumption, allowing varying weights for different closures and separating closure effects for land and ocean. For the surface, CWRF includes a comprehensive multilevel turbulence upper ocean model (UOM) to resolve the transit air-sea interactions (Ling et al., 2011, 2015). The UOM has recently been modified to improve shallow lake modeling (L. Sun et al., 2020).

Most relevant to this study is that CWRF couples the state-of-the-art CSSP to represent land surface interactions with the atmosphere. For consistent coupling, CWRF incorporates an advanced dynamic-statistical parameterization of land surface albedo (Liang, Xu, et al., 2005) together with realistic surface characteristics distributions (Liang, Choi, et al., 2005; M. Xu et al., 2014). In contrast to most land surface models which typically capture soil-moisture transport only in the vertical direction, CSSP is built with a 3-D conjunctive surface-subsurface flow model to represent the lateral and subgrid soil-moisture transport, the surface flow routing and interaction with subsurface flow, and the topographically controlled baseflow (Choi et al., 2007, 2013; Choi & Liang, 2010; Yuan & Liang, 2011a). It has been demonstrated with advanced skills in predicting detailed terrestrial hydrology variations, soil temperature and moisture distributions, and land-atmosphere interactions, as well as expanded applications (Gan et al., 2015; Ji et al., 2017; Yuan et al., 2018). Figure 1 shows the 30-km flow directions resolved by CSSP in the Great Lakes watershed. They catch runoffs throughout the watershed and route the accumulated stream flows into the Great Lakes. In particular, a mass conservative method was used to treat river discharge from CSSP at 30-km grids to FVCOM unconstructed nodes at ~1 km near shore. For a given lake, all land grids within its entire watershed or basin contribute runoff to the lake. The CSSP has already represented at each grid the accumulated runoff from all contributing upstream grids within the basin (Choi et al., 2013). Such accumulative runoff is routed as one representative river discharge to the nearest and deepest FVCOM node within each CWRF coastal grid. To avoid double counting, if a surrounding-lake grid is located upstream of another coastal grid, it is discarded.

2.2 Lake and Ice Joint Component (FVCOM)

FVCOM is a free-surface, unstructured-grid, finite-volume, 3-D primitive equation coastal ocean model originally developed and continuously improved by C. Chen et al. (2003, 2006). The unstructured triangular mesh used in horizontal gridding is highly advantageous for adequately representing the complex geometry such as estuaries and coastlines. The primitive equations are discretized vertically using the terrain-following sigma coordinate. FVCOM is integrated using an explicit mode-splitting approach, which solves the barotropic and baroclinic modes separately with a short and a longer time step.

Recently, FVCOM has incorporated an unstructured-grid, finite-volume sea ice model (UG-CICE) developed by Gao et al. (2011) based on Los Alamos Community Ice CodE (CICE, Hunke & Lipscomb, 2008). It predicts spatiotemporal variations of ice thickness by combining a thermodynamic model to compute ice local growth rates, an ice dynamics model to describe ice pack velocity, a transport model to determine the advection of state variables, and a ridging parameterization to represent the mechanical processes. Solving the lake ice dynamics is computationally intensive, approximately tripled the FVCOM computing time. Since the ice cover over the Great Lakes is commonly confined within the closed basins, only the ice thermodynamic part (Bitz & Lipscomb, 1999) is activated in this study to calculate the ice growth rate with a remapping scheme (Lipscomb, 2001) to account for ice transport in thickness space. Adding FVCOM increased the overall CWRF computing time by ~40%.

2.3 Coupler (CPL7) and CWRF-FVCOM Coupling

A key issue for constructing a coupled system is how to connect component models of distinct infrastructures like CWRF and FVCOM. The traditional approach is to treat these components as internal routines that can directly call each other. This, however, lacks flexibility and extensibility and makes it difficult to incorporate other components. The modern approach for building a modular and flexible coupled system is to use advanced coupling techniques (i.e., couplers) that interact and coordinate separate components by facilitating mesh decomposition, data transfer, grid mapping, time scheduling, and other functionalities (Collins et al., 2005; Craig et al., 2005, 2012; Jacob et al., 2005; Larson et al., 2005; Valcke, 2006). CPL7 is a popular coupler (Craig et al., 2012), developed by the National Center for Atmospheric Research. For example, it was used to nest WRF within the Community Earth System Model for online downscaling from coarse global circulation to refined cyclogenesis over the U.S. Southern Great Plains (He et al., 2013). It was also adopted for the Regional Arctic System Model to link five components in representing interactive processes in the Arctic (Cassano et al., 2017; Hamman et al., 2016).

We choose CPL7 to couple CWRF with FVCOM, which has already integrated with UG-CICE (Gao et al., 2011) and other models for future expansion (Figure 2a). Each component has predefined processor groups and distinct physical models that can run sequentially or concurrently and exchange data only with the coupler. CPL7 not only contains a top driver that runs on all processors and controls component sequencing, processor decomposition, and data communication (Craig, 2014; Craig et al., 2012) but also provides key functionalities such as data mapping (i.e., interpolated in same processors), rearranging (i.e., data transfer among different processor groups), and merging using the Model Coupling Toolkit (MCT) library.

The coupled system integrates sequentially (Figure 2b), that is, FVCOM and CWRF take turns integrating and import (export) data are rearranged from (to) CPL7 before (after) each corresponding coupling time interval. The exchange data and mesh decompositions (CWRF: rectangular; FVCOM: unstructured triangular) are both represented in MCT data types. Mapping between these two distinct meshes is achieved through sparse matrix interpolation with the weights precomputed using the Spherical Coordinate Remapping Interpolation Package (SCRIP, Jones, 1998).

Moreover, standard coupling interfaces are particularly designed following the Earth System Modeling Framework structure to connect the top-level driver with the model libraries which are compiled from modified CWRF and FVCOM codes. The major modifications include (1) replacing the default root Message Passing Interface communicator with the one directly allocated from the driver, (2) synchronization of time managers between models and the driver, (3) construction of derived data types storing variables from the coupler, and (4) update of required forcings or surface boundary conditions directly from new data structures.

3 Experiment Design and Evaluation Data 3.1 Experiment Design

Table 1 summarizes the primary CWRF-FVCOM configuration in this study. The CWRF computation domain (Figure 3) centered at (37.5°N, 95.5°W) covers the whole contiguous U.S. and adjacent regions under the Lambert conformal map projection. The model consists of 197 (west-east) × 139 (south-north) horizontal grids at 30-km spacing and 36 vertical layers up to 50 hPa. Along the four edges of the domain, a width of 14 grids is identified as the buffer zones where varying lateral boundary conditions (LBCs) are specified using a dynamic relaxation technique (Liang et al., 2001). This domain design has demonstrated superior downscaling skill for U.S. climate (L. Chen et al., 2016; Liang et al., 2001, 2007, 2012; Liu et al., 2015; Yuan & Liang, 2011b). The CWRF control physics configuration includes radiation—GSFCLXZ (Chou et al., 2001; Chou & Suarez, 1999); land surface—CSSP (Choi et al., 2007); planetary boundary layer—CAM (Holtslag & Boville, 1993); cumulus—ECP penetrative convection (Qiao & Liang, 2015, 2016, 2017) plus UW shallow convection (Park & Bretherton, 2009); microphysics—GSFCGCE (Tao et al., 2003); ocean—UOM (Ling et al., 2015); and Lake (excluding Great Lakes grids)—LISSS (Subin et al., 2012). For more details refer to Liang et al. (2012).

Table 1 Summary of CWRF-FVCOM Configuration

CWRF
Domain			Contiguous United States and adjacent regions, centered at (37.5°N, 95.5°W)
			Horizontal resolution: 30 km (196 × 139 grids)
			Vertical resolution: 36 levels, top at 50 hPa
			Buffer zone width: 420 km (14 grids)
Physics	Cloud-aerosol-radiation		GSFC (Chou et al., 2001; Chou & Suarez, 1999; Liang & Zhang, 2013)
	Land surface		CSSP (Choi et al., 2007, 2013; Choi & Liang, 2010; Dai et al., 2003, 2004)
	Planetary boundary layer		CAM (Holtslag & Boville, 1993)
	Cumulus		Deep convection: ECP (Qiao & Liang, 2015, 2016, 2017)
			Shallow convection: UW (Park & Bretherton, 2009)
	Microphysics		GSFCGCE (Tao et al., 2003)
	Ocean		UOM (Ling et al., 2015)
	Lake (default)		LISSS (Subin et al., 2012)
FVCOM
Domain			Great Lakes (Superior, Michigan, Huron, Erie, and Ontario)
			Horizontal resolution: 1–8 km; (nodes:14481; elements: 27194)
			Vertical resolution: 31 layers
Physics	Turbulent closure		Horizon: Smagorinsky (Smagorinsky, 1963)
			Vertical: MY-2.5 (Mellor & Yamada, 1982)
	Surface flux scheme		Open water (FVCOM): COARE 2.6 (Fairall et al., 1996)
			Ice surface (UG-CICE): J99 (Jordan et al., 1999)
	Ice	Thermodynamics	Bitz and Lipscomb model (Bitz & Lipscomb, 1999)
		Dynamics	Turn off
		Thickness distribution	Remapping scheme (Lipscomb, 2001)
Channel			K	YM	A	B	Zc	Wt
	St. Marys		824.7	181.43	1.3	0	0	1.0
	St. Clair		56.4	165.83	2.01	0.5	−0.02	0.5
	Detroit		70.7	165.94	2.1	0.5	−0.01	1.0
	Niagara		616.7	169.75	1.67	0	0	1.0
	St. Lawrence		555.82	69.50	1.65	0	0	1.0
CPL7
Coupling frequency			1 hr
Exchange Data	CWRF to CPL7		Meteorological fields at the lowest model level (wind, pressure, temperature, and relative humidity)
			Radiations (downward shortwave and longwave radiative fluxes)
			Precipitation, evaporation, and runoff
	FVCOM to CPL7		Lake surface temperature (LSTa) and lake ice cover (LIC)

aLST = WST × (1 − LIC)+IST × LIC, where WST (IST) is water (ice) surface temperature.

The FVCOM computation domain (Figure 1) covers all Great Lakes. It consists of 14,481 unstructured nodes and 27,194 triangular elements at varying horizontal spacing (1–8 km) and 31 vertical terrain-following sigma layers. Bathymetry data provided from the National Geophysical Data Center are interpolated to the model grids. The horizontal and vertical mixing are parameterized respectively using Smagorinsky (1963) and Mellor and Yamada (1982) level 2.5 turbulence closure scheme. Surface sensible and latent heat fluxes are parameterized in open water periods by the FVCOM's own COARE2.6 algorithm (Fairall et al., 1996), which agrees with eddy covariance measurements better than other schemes (Charusombat et al., 2018), and during ice-covered periods by the J99 scheme built-in UG-CICE (Jordan et al., 1999). Solar radiation penetrates into the upper water layers following a two-band scheme (Kraus, 1972).

The integration time intervals differ largely among components and their own processes. CWRF at 30-km grid spacing uses 120 s for the dynamics, microphysics, and convection, 10 min for the land/ocean surface and planetary boundary layer, and 30 min for the radiative transfer. In contrast, FVCOM at finer unstructured grids (1–8 km) uses 5 and 30 s as the external and internal time steps to, respectively, integrate the 2-D and 3-D processes. The CWRF-FVCOM allows varying specifications of the coupling frequencies among the different components, and this study uses 1 hr for each component exchanging data with CPL7. FVCOM is configured to receive inputs of lowest-level meteorological fields (pressure, temperature, humidity, and wind), surface downwelling radiation fluxes, precipitation, evaporation, and runoff from the coupler. In turn, the FVCOM predicted lake surface temperature (LST) and lake ice cover (LIC) are sent back to the coupler as the surface boundary conditions for CWRF. In a lake grid with ice cover, LST = WST × (1 − LIC)+IST × LIC, where WST is the open-water surface temperature predicted by FVCOM and IST is the temperature at the ice/air interface predicted by UG-CICE. As such, CWRF and FVCOM/UG-CICE calculate separately surface sensible and latent heat fluxes and wind stresses, which are used independently for the vertical transfer processes into the atmosphere and the water.

This above coupling strategy has been adopted in the design of several coupled modeling systems (Sitz et al., 2017; Xue et al., 2017). An important exception is our use of precipitation, evaporation, and runoff as additional exchange fields to account for the water balance. Such surface-atmospheric coupling based on a dynamic representation of meteorological conditions shows advantages in reducing surface forcing uncertainty in hydrodynamic modeling (Xue et al., 2015). In an off-line test, FVCOM produces a more realistic surface temperature when using the net heat flux calculated by its built-in bulk formula than that provided from CWRF. A similar result has been reported (Xue et al., 2015). Note that FVCOM receives evaporation from the coupler as the external forcing together with precipitation and lateral river discharge. Doing so not only minimizes structural changes to FVCOM but also ensures the water balance between surface and atmosphere exchanges. Moreover, additional online tests, in which CWRF receives FVCOM computed heat fluxes, presents minor differences as compared to the current coupling. In the future, we may adopt the FVCOM bulk formula into CWRF for a consistent surface heat exchange between the two.

We integrated CWRF-FVCOM continuously from 1 March 1998 to 31 December 2015. The CWRF initial and LBCs were generated from the European Centre for Medium-Range Weather Forecasts Interim reanalysis (ERI, Dee et al., 2011). Sea surface temperature was dynamically predicted using UOM (Ling et al., 2015) with daily nudging of the observational analysis (Banzon et al., 2016; Reynolds et al., 2007). FVCOM was initialized in a motionless state with a uniform water temperature of 2°C throughout the Great Lakes. The first 9 months were regarded as spin-up (adequate for dimictic lakes) and not used for subsequent analysis. In parallel, a baseline integration (CWRF-LISSS) was conducted using the CWRF control configuration. This run was identical to CWRF-FVCOM except that LISSS (Gu et al., 2015; Subin et al., 2012) replaced FVCOM in the Great Lakes.

3.2 Evaluation Data

Table 2 summarizes the key observational data sets and their sources used in this study as the reference for model evaluation. The Great Lakes Surface Environmental Analysis (Schwab et al., 1992) provides daily digital WST maps derived from National Oceanic and Atmospheric Administration (NOAA) polar-orbiting satellite imageries from 1995 onward. The National Data Buoy Center operates in situ routine observations at nine buoys (Figure 4), including water temperature at 0.5-m depth (treated here as WST) and surface air temperature. The National Ice Center produces an ice cover analysis from 1998 onward at intervals of 3–7 days, which are linearly interpolated into daily maps over each of the five Great Lakes.

Table 2 Summary of the Evaluation Data Set Used in This Study

Data set	Variable	Evaluation period	Time interval	Grid resolution	Reference	Source
Great Lakes Surface Environmental Analysis	WST	1999–2015	Daily	0.01°	Schwab et al. (1992)	https://coastwatch.glerl.noaa.gov/thredds/Satellite/glsea/glsea_catalog.html
National Data Buoy Center Station Data	WST; surface air temperature	2002–1,015	Daily	—	—	https://www.ndbc.noaa.gov/
Great Lakes Ice Analysis Products	Ice cover	1999–2015	Daily	0.01–0.03°	—	https://www.natice.noaa.gov/products/great_lakes.html
Modern-Era Retrospective analysis for Research and Applications, Version 2	Surface air temperature	1999–2015	Monthly	0.625° × 0.5°	Gelaro et al. (2017)	https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/data_access/
Climatology-Calibrated Precipitation Analysis	Precipitation	2002–2015	Monthly	0.125° × 0.125°	Hou et al. (2014)	https://para.nomads.ncep.noaa.gov/pub/data/nccf/com/ccpa/para/
Great Lakes Environmental Research Laboratory Monthly Hydrometeorological Database	NBS components	1999–2015	Monthly	—	Hunter et al. (2015)	https://www.glerl.noaa.gov/pubs/tech_reports/glerl-083/UpdatedFiles/
Coordinated monthly mean Lakewide average water levels	Water level	1999–2015	Monthly	—	—	http://lre-wm.usace.army.mil/ForecastData/GLHYD_data_metric.csv

For surface air temperature, the reference comes from the latest Modern-Era Retrospective analysis for Research and Applications Version 2 at ~50-km grid produced by the National Aeronautics and Space Administration (Gelaro et al., 2017). For precipitation, the reference is from the NOAA Climatology-Calibrated Precipitation Analysis merging daily gauge and 6-hourly Stage IV data at 5-km grid (Hou et al., 2014). For the NBS components, the reference is derived from monthly estimates of overlake evaporation, overland and overlake precipitation, and runoff included in the NOAA Great Lakes Environmental Research Laboratory monthly hydrometeorological database (Hunter et al., 2015). For the lake water level, the reference is based on coordinated monthly data from the U.S. Army Corps of Engineers.

4 Results 4.1 Lake Conditions

Given that observed IST data are not available, below we evaluate WST rather than LST. Figure 5 compares the 1999–2015 mean annual cycles of daily WST averaged over each of the five Great Lakes between CWRF-FVCOM and CWRF-LISSS, including also their mean biases, root-mean-square errors, and interannual correlations relative to observations. The WST annual cycle can be divided to a warming phase from March to mid-August and a cooling phase for the rest of the year. Deeper lakes typically exhibit cooler water temperatures in summertime. CWRF-FVCOM accurately capture the annual cycle in both phase and magnitude, with overall small errors of 0.66°C, 0.59°C, 0.82°C, 0.99°C, and 0.73°C for Lakes Superior, Michigan, Huron, Erie, and Ontario, respectively. In contrast, CWRF-LISSS produces much larger errors (3.1°C, 1.97°C, 2.07°C, 2.38°C, and 1.75°C) for all five lakes. It simulates too rapid spring warming that leads to significantly warmer summer (especially in Lake Superior by up to 6°C) and autumn (except for Lake Erie). The result indicates the inefficiency of 1-D lake models that tend to predict earlier stratification at large and deep lakes due to lack of horizontal processes (Bennington et al., 2014; Martynov et al., 2012; Xiao et al., 2016). In particular, both horizontal and vertical velocity shears determine lake mixing processes. LISSS simply treats vertical mixing as a function of surface wind with a factor depending only on lake depth. This is inadequate to resolve spatial variability in lake mixing processes. Moreover, horizontal heat variations in LISSS cannot be redistributed due to the lack of advective transport through wind- or density-driven circulations. On the other hand, FVCOM with a 3-D hydrodynamic representation significantly improves the simulation of the thermal structure in the Great Lakes characterized with a wide range of water depths.

Figure 5 also shows that CWRF-FVCOM systematically improves over CWRF-LISSS, with much smaller mean biases and overall errors as well as higher interannual correlations in all lakes. Figure S1 in the supporting information further illustrates that CWRF-FVCOM captures realistic lake-wide average WST interannual variations. For example, it well reproduces the lower WST values in 2003, 2009 and 2014, each characterized by larger ice cover in winter, as well as the stronger variability during 2012–2014. Figure S2 compares CWRF simulations against available measurements at the nine buoys, which locate across the Great Lakes with varying depths. The result again reveals that CWRF-FVCOM outperforms CWRF-LISSS in simulating WST daily to interannual variations during 1999–2015 both at shallow and deep lake buoys.

The Great Lakes show strong WST spatial variations across the wide horizontal scales and complex bathymetries. Figure 6 compares CWRF 1999–2015 averaged seasonal mean WST distributions against observations. CWRF-FVCOM shows excellent skills in simulating both the seasonal and spatial variations with detailed characteristics. It well reproduces the cold WST (<4°C) in winter, quick warming of nearshore waters in spring, distinct temperature differences between nearshore and offshore waters in summer, and coast-to-offshore cooling in autumn. In contrast, CWRF-LISSS is less skillful, especially in spring and summer. The rapid WST warming in spring makes cold waters to reach 4°C much more quickly, which leads to an earlier strong stratification and consequently causes a significant warm bias in summer. It also overestimates WST in autumn except for Lake Erie.

The effect of ice cover on lake thermal and regional climate has been well recognized (Austin & Colman, 2007; Notaro, Bennington, & Vavrus, 2015; Vavrus et al., 2013). For example, extended ice cover can significantly reduce winter lake-effect precipitation and decline summer water temperatures due to the retarded stratification. Figure 7 compares CWRF LIC annual cycle and interannual variations against observations over each of the Great Lakes. CWRF-FVCOM realistically captures the seasonal to interannual variations of each lake, as well as the contrasts among the lakes, for example, more extensive ice cover in Lake Erie than Ontario. It also produces accurate ice onset in early December and slightly stronger ice melting in March. Its simulated extensive cover in severe winters (2003, 2009, and 2014) and limited cover in mild winters (2006 and 2009) agree well with observations. On the other hand, CWRF-LISSS, which defines LIC as the mass fraction of lake water frozen in the top layer, substantially underestimates the ice cover magnitude of both the annual cycle and interannual variability in all lakes. Notably, CWRF-FVCOM also significantly improves over WRF-LISSS in simulating the ice cover variations during 2012–2014 (Xiao et al., 2016). However, it is beyond the scope of this study to determine which factors may contribute to LISSS biases or improve its ice process representation.

Figure 8 compares CWRF-FVCOM and CWRF-LISSS simulated vertical temperature profiles at Station 45006 in Lake Superior. In observations (Austin & Colman, 2008), the annual length of positively stratified season (i.e., WST > 3.98°C) is around 170 days with the onset in mid-June. As coupled with CWRF, FVCOM faithfully reproduces the thermal features, while LISSS produces an early onset and overextended length of the stratification. The upper hypolimnion depth simulated by FVCOM is between 30 and 40 m during summertime, which agrees well with the estimate of W. Xu et al. (2019) based on in situ measurements. However, LISSS produces an insufficiently sharp temperature gradient in the metalimnion layer, indicating its inefficient vertical mixing. It simulates too rapid spring warming due to the combined effects of its lacking horizontal processes and underestimating ice cover (Figure 7). This results in earlier ice break up and hence earlier start of stratification, causing the water to be overheated. Therefore, LISSS predicts a longer period of stratification, a deeper thermocline, and later autumn cooling. FVCOM also outperforms LISSS in representing interannual variability of the vertical thermal regimes under the joint effects of climate variation and positive ice albedo feedback. For example, it realistically reproduces the delayed (advanced) and weak (strong) stratification in cold (warm) years such as 2009 and 2014 (2006 and 2012), leading to cooler (warmer) water temperature during summertime. Figure S3 compares the simulated monthly mean vertical profiles at Station 45007 in Lake Michigan in 2011, indicating again that CWRF-FVCOM produces the stratification onset and peak closer (than CWRF-LISSS) to observations (Xiao et al., 2016).

The long-term currents in the Great Lakes are mainly regulated by wind stress (wind-driven) and surface heat flux (density-driven) in the context of complex bathymetries (Beletsky et al., 1999; Pickett, 1980). In particular, the wind-driven circulation, controlled by wind stress and horizontal pressure gradient, dominates in the cold season, whereas the density-driven circulation, regulated by horizontal pressure gradient and Coriolis force, becomes comparable in the warm season. Beletsky et al. (1999) and Bai et al. (2013) gave a comprehensive overview of seasonal climatological circulation patterns based on observations and simulations. Overall, the mean depth-averaged currents produced by CWRF-FVCOM (Figure 9) share similar patterns to the climatological circulation map in Bai et al. (2013). In winter, CWRF-FVCOM reasonably reproduces the general cyclonic circulations in Lakes Superior, Michigan, and Huron, which are mainly driven by the prevailing overlake westward or northwestward winds. It also realistically captures the strong currents such as in the east of the Keweenaw Peninsula and along the east coast of Lake Michigan. In summer, the currents become weaker and more complex due to the decrease of wind speeds and enhanced baroclinic effects respectively. Lake Michigan is characterized by a (an) cyclonic (anticyclonic) circulation in northern (southern) basin with deeper (shallower) depth. An obvious anticyclonic (cyclonic) gyre dominates most parts on Lake Erie (Ontario). In general, our model simulates similar seasonal circulation patterns to observations (Beletsky et al., 1999) and other model studies (Bai et al., 2013), indicating the credibility of CWRF-FVCOM. Some differences in simulations exist due to discrepancies such as model configuration, spatial resolution, and time period.

4.2 Regional Climate Effects

Figure 10 compares CWRF-FVCOM and CWRF-LISSS geographic distributions of 1999–2015 averaged seasonal mean surface air temperature biases over the Great Lakes region, while Figure 11 depicts their monthly mean differences averaged separately over the lake, land, and basin grids (see Figure 1 for the boundaries). In winter, CWRF-LISSS produces a cold bias over Lake Superior, which CWRF-FVCOM reduces by ~1°C. For the other seasons, CWRF-LISSS strongly overestimates overlake air temperature, which CWRF-FVCOM significantly improves by up to 4.7°C, 3.5°C, 2.9°C, 2.3°C, and 2.8°C over Lakes Superior, Michigan, Huron, Erie, and Ontario, respectively. The reduction is mainly attributed to the improved WST simulation, which substantially decreases the turbulent heat fluxes (Figures S4 and S5), especially the latent heat (by up to 45, 42, 34, 61, and 36 W m⁻² over the respective lakes). The systematic CWRF overestimates of air temperature over land (Liang et al., 2012) mainly near the lake shoreline are also reduced, especially in June by ~0.5°C. Over land, skin surface temperature increases more than the overlying air temperature in warm seasons (Figure S6), enhancing local sensible heat flux (Figure S5). Moreover, low-level winds show a geostrophic response to the air temperature changes, with an obvious anticyclonic structure in summer and autumn due to cooling (Figure 10). Figure S7 further compares simulations against in situ measurements (~4 m above lake surface) at the nine buoys. Clearly, CWRF-FVCOM produces more realistic air temperature variations than CWRF-LISSS, especially during the stratification period at deep lake locations.

Figure 12 compares the geographic distributions of 1999–2015 averaged monthly mean CWRF-FVCOM minus CWRF-LISSS differences in surface latent and sensible heat fluxes as well as precipitation from December to March over the Great Lakes region, while Figure 13 depicts the corresponding precipitation annual cycle differences averaged separately over the lake, land, and basin grids. The overlake responses are overall weaker in precipitation than air temperature. Significant precipitation changes occur only over Lake Eire in July. During late autumn and winter, relatively cold, dry air masses crossing over the warm water get heated and moistened, resulting in lake-effect precipitation on the downwind shores of the lakes. Both CWRF-FVCOM and CWRF-LISSS can reproduce the general patterns of precipitation, such as its southeast to northwest gradient and distinct upwind to downwind contrast across the lakes (Figure S8). CWRF tends to overestimate precipitation on the downwind shores, especially for Lakes Huron, Erie, and Ontario. Similar overestimation was identified in WRF downscaling (Gula & Peltier, 2012). However, the reference data contain large uncertainty in the Great Lakes region (Gula & Peltier, 2012; Xiao et al., 2016), and we cannot separate model deficiencies from data errors.

Nonetheless, the CWRF-FVCOM minus CWRF-LISSS differences in cold season precipitation over the lakes and downwind regions can be interpreted mainly from the thermodynamic perspective (Figure 12). Turbulent heat fluxes are usually upward in cold seasons. In December, CWRF-FVCOM relative to CWRF-LISSS simulates smaller overlake turbulent fluxes except for Lake Erie's central basin and several nearshore lake grids, which causes overlake and downwind precipitation declines up to 0.4 mm/day. When getting colder, the flux difference reverses its sign in Lake Superior, leading to increased local precipitation (~0.1 mm/day). CWRF-LISSS produces colder air temperature than CWRF-FVCOM over Lake Superior in January and February (Figure 11) in spite of its substantial underestimation of ice cover (Figure 7). As the stratification occurs, intensifies, and collapses (April to November), CWRF-FVCOM generally reduces precipitation from CWRF-LISSS over each lake basin. The reduction peaks mostly in June or July, ranging in 0.25–1.0 mm/day over the lake and in 0.19–0.67 mm/day over land of the basin. The surface and air temperature difference can measure atmospheric stability, favoring a stable or unstable condition when it is negative or positive (Holman et al., 2012). Therefore, the precipitation reduction by CWRF-FVCOM from CWRF-LISSS is mainly due to decreased moisture and heat fluxes as well as enhanced atmospheric stability (Figure S6) over the lake.

4.3 Lake Water Levels

Hydrodynamic models are used commonly to predict lake water levels for short-term fluctuations caused by varying meteorological conditions but are rarely to do so for long-term variations driven primarily by NBS and connecting channel flows that are difficult to determine. In both cases, hydrodynamic models run off-line as driven by forcing inputs from various data sources, but they are decoupled from regional climate models that provide hydrometeorological forcings. To bridge this gap, our fully coupled system transfers CWRF simulated precipitation, evaporation, and runoff to FVCOM for surface water budgets and lateral inflows. The system accounts for connecting channel flows using the stage-fall-discharge relationship adopted from the Coordinated Great Lakes Regulation and Routing Model (Clites & Lee, 1998; Quinn, 1978) as [Image Omitted. See PDF]where Z_u and Z_d are predicted lake-wide mean water levels of upstream and downstream lakes. K, YM, A, B, and Z_c are the equation coefficient, channel invert, depth exponent, fall exponent, and fall adjustment constant, respectively. W_t is the weighting factor of upstream lake. Table 1 lists the values of these parameters used in this study. All elevations Z_x are referenced to the International Great Lakes Datum (IGLD-85) in the unit of meter. Lakes Michigan and Huron (including Georgian Bay) are considered as one lake since they are hydraulically connected. Diversions and ice retardation are included using the monthly climatology suggested in the above model. Other minor components such as groundwater and consumptive use are neglected.

Figure 14 compares CWRF-FVCOM simulated with the reference's monthly mean NBS and water level variations during 1999–2015. Note that the reference estimates contain large uncertainties, especially for NBS (DeMarchi et al., 2009). Overall, the coupled system can reasonably reproduce the NBS variability. It generally captures the high NBS values in spring to early summer and lower values for the rest of the year (Neff & Nicholas, 2005) as well as strong interannual variations. In particular, the simulated water level of Lake Superior is quite realistic before 2011, followed by underestimation. The water level of Lakes Michigan-Huron is overestimated in 2002–2006 and underestimated in 2011–2015. The simulation tends to be worse for Lakes Erie and Ontario, having smaller correlations (0.58 and 0.59) as compared with Lakes Superior and Michigan-Huron (0.85 and 0.67). As discussed below, water level prediction critically depends the accuracy of NBS variations in each lake. However, it is still a challenge to improve NBS simulation due to the difficulty in accurately representing its contributing components, especially evaporation and runoff that even lack good quality observations as the reference (Notaro, Bennington, & Lofgren, 2015).

Figure 15 compares CWRF-FVCOM simulated with the reference's 1999–2015 averaged annual cycles of NBS and its components (precipitation, evaporation, and runoff) as well as water level for each lake. The model generally represents seasonal variations of lake water level as it overall well captures the annual cycle of NBS. The annual cycle simulation is quite realistic for both water level and NBS of Lakes Superior and Michigan-Huron, despite the deficiencies noted earlier (e.g., the cancelation between overestimation and underestimation for the latter two lakes). In these lakes, precipitation prediction is reasonably accurate, considering the observational uncertainty range. Model disparities from the reference are larger for overlake evaporation and lateral runoff, both of which are associated with large data uncertainties. Model to reference departures are more obvious in Lakes Erie and Ontario for NBS and all components, probably because of their smaller lake sizes implied with larger data uncertainties. Lake Erie is identified with notably larger seasonal variations in both evaporation and runoff than Lakes Superior and Michigan-Huron and greater simulated NBS departures from the reference. In contrast, Lake Ontario has the largest runoff contribution among all five lakes. Inadequate representation of its outflow is partly responsible for limiting model prediction skill of this lake's water level (Dupont et al., 2012). In fact, the channel flows through the St. Marys and St. Lawrence Rivers are regulated by International Joint Commission depending on the real-time water level. This regulation has not been incorporated into the model. Furthermore, our off-line tests show that the interannual trend of unregulated lake water level is highly sensitive to the specification of depth exponent A. Therefore, there is still large room for improvement to water level prediction by more accurate representation of NBS' components along with optimal specification of connecting channel flows' parameters.

5 Conclusion and Discussion

There is a great need for an integrated modeling system to improve the understanding and prediction of land-lake-atmosphere interactions and water balance over the Great Lakes region (Sharma et al., 2018). Although 3-D hydrodynamic models significantly improve over 1-D models in representing lake thermal structures, it is challenging to couple them with regional climate models because of technical difficulties in efficiently linking the two components that differ dramatically in design infrastructures. Moreover, hydrodynamic models are rarely implemented to predict the long-term water level variations of the Great Lakes because of the necessity to account for NBS and connecting channel flows that are difficult to determine accurately. This study has developed an interactive lake-atmosphere-hydrology modeling system by coupling regular-grid CWRF with unstructured-grid FVCOM using the CPL7 coupler to predict atmosphere-watershed interactions over the Great Lakes region. The coupled CWRF-FVCOM is compared with CWRF-LISSS for historical simulations during 1999–2015 against observations to investigate how 3-D hydrodynamic processes improve predictions of lake thermal structures and regional climate patterns and also to explore its ability to predict seasonal-interannual water level variations.

CWRF-FVCOM outperforms CWRF-LISSS in simulating water surface temperature, ice cover, and vertical thermal structure at seasonal to interannual scales for all five Great Lakes. It also realistically reproduces the seasonal lake circulation patterns. The improved lake conditions in turn significantly reduce air temperature overestimates in warm seasons by up to 4.7°C, 3.5°C, 2.9°C, 2.3°C, and 2.8°C over Lakes Superior, Michigan, Huron, Erie, and Ontario respectively, corresponding mainly to decreases in surface latent heat flux by up to 45, 42, 34, 61, and 36 W m⁻². The low-level wind shows a geostrophic response to the temperature reduction, with an obvious anticyclonic structure in summer and autumn over the Great Lakes region. Meanwhile, precipitation is generally reduced, although not as significantly as air temperature. Significant precipitation reductions occur only over Lake Eire in July. In cold months, precipitation changes over the lakes and downwind are highly related to changes in surface turbulent heat fluxes. In warm seasons, the general reduction of precipitation over the lake basin is mainly related to the decreased moisture and heat fluxes as well as the enhanced atmospheric stability over the lake. While it is difficult to make a direct comparison with other off-line (Bai et al., 2013) or online (Xue et al., 2017) modeling studies, our results show at least equally, if not more, skillful performance against available measurements. The key to improve modeling the Great Lakes (especially for Lake Superior) is to incorporate a 3-D hydrodynamic model to account for the horizontal processes unresolved in a 1-D lake model. As a result, almost all 3-D modeling efforts over one or all of the Great Lakes can achieve similarly good performance.

For the water balance over the Great Lakes basin, CWRF is able to generally capture NBS seasonal to interannual variations. The coupled CWRF-FVCOM, when adopting a stage-fall-discharge equation to account for connecting channel flows, reasonably represents long-term water level variations, especially for Lakes Superior and Michigan-Huron. However, refinements are needed to improve modeling NBS components, especially evaporation and runoff, as well as channel flows, especially for Lakes Erie and Ontario. It is also helpful to test the importance of ice dynamics omitted in this study. Consequently, the coupled system can be operated to predict lake thermal structure and water level variations not only in response to current climate anomalies but also as regional consequences from global climate change.

Considering that our main focus was on developing a coupled system that can be applied on climate scales, the CWRF resolution will have to be limited (30 km in this study, albeit currently testing at 10 km) to meet the computational demand. This was also one of the reasons why we sought the unconstructed grid FVCOM to resolve the Great Lake dynamical processes. Ideally, the CWRF's hydrologic component (CSSP) would be run at a fine resolution similar to that of FVCOM near shore (~1 km), thus making natural connections from the routing rivers to the respective lake nodes. This concept was in our original design for future implementation of the next generation regional Earth system model (5–10-km atmosphere, 1-km land, and 1- to 10-km lakes and oceans). On the other hand, the realistic simulations of water surface temperature, vertical thermal structure, and ice cover over the Great Lakes as well as seasonal circulation patterns over the larger region are in favor of the effectiveness of our current mesh configuration and CWRF-FVCOM coupling. Given the promising performance over the Great Lakes, ongoing efforts are to refine CWRF grid resolution and physics representation (especially for hydrology), build the unstructured meshes for the U.S. coastal oceans, and develop subseasonal-to-seasonal regional climate prediction infrastructure. Ultimately, the coupled CWRF-FVCOM system can be applied to predicting water levels in the Great Lakes and U.S. coastal oceans as well as their interactions with regional climate and watersheds.

Acknowledgments

The research was supported by U.S. National Science Foundation Innovations at the Nexus of Food, Energy and Water Systems under Grants EAR1639327 and EAR1903249. The simulations and analyses were conducted on supercomputers, including the Computational and Information Systems Lab of the National Center for Atmospheric Research and the Maryland Advanced Research Computing Center's Bluecrab. We thank Chao Sun for helping with the CWRF model configuration.

Data Availability Statement

Observational data and model simulations used in this study are accessible online (at ftp://earthserver.umd.edu/publish).