Insights on Simulating Summer Warming of the Great Lakes: Understanding the Behavior of a Newly Developed Coupled Lake‐Atmosphere Modeling System

The Laurentian Great Lakes are the world's largest freshwater system and regulate the climate of the Great Lakes region, which has been increasingly experiencing climatic, hydrological, and ecological changes. An accurate mechanistic representation of the Great Lakes thermal structure in Regional Climate Models (RCMs) is paramount to studying the climate of this region. Currently, RCMs have primarily represented the Great Lakes through coupled one‐dimensional (1D) column lake models; this approach works well for small inland lakes but is unable to resolve the realistic hydrodynamics of the Great Lakes and leads to inaccurate representations of lake surface temperature (LST) that influence regional climate and weather patterns. This work overcomes this limitation by developing a fully two‐way coupled modeling system using the Weather Research and Forecasting model and a three‐dimensional (3D) hydrodynamic model. The coupled model system resolves the interactive physical processes between the atmosphere, lake, and surrounding watersheds; and validated against a range of observational data. The model is then used to investigate the potential impacts of lake‐atmosphere coupling on the simulated summer LST of Lake Superior. By evaluating the difference between our two‐way coupled modeling system and our observation‐driven modeling system, we find that coupled‐lake atmosphere dynamics can lead to a higher LST during June‐September through higher net surface heat flux entering the lake in June and July and a lower net surface heat flux entering the lake in August and September. The unstratified water in June distributes the entering surface heat flux throughout the water column leading to a minor LST increase, while the stratified waters of July create a conducive thermal structure for the water surface to warm rapidly under the higher incoming surface heat flux. This research provides insight into the coupled modeling system behavior, which is critical for enhancing our predictive understanding of the Great Lakes climate system.

the Finite Volume Community Ocean Model (FVCOM) (Chen et al., 2006(Chen et al., , 2013. Sun et al. (2020) followed up by two-way coupling the Climate-Weather Research and Forecasting to an FVCOM-based 3D lake and ice model. Xue et al. (2022) and Kayastha et al. (2022) later built upon the framework from Xue et al. (2017) and derived projections for the Great Lakes under climate change. However, as far as we know, Xue et al. (2017) and Sun et al. (2020) are still the only two studies dedicated to documenting the two-way coupling of an RCM to a 3D lake model for the Great Lakes.
Furthermore, due to the rarity of two-way coupled RCMs with 3D lake models for the Great Lakes, a complete mechanistic understanding of the impact of lake-atmosphere coupling on the simulation of Great Lakes LST is lacking. Unlike in historical simulations where available observational data can be used to supplement models or even circumvent model limitations, future projections for the Great Lakes LST rely completely on the results from RCM and lake model including the lake-atmosphere coupling between them. The role of such an interactive process between the atmosphere and lake model on the simulated LST has not been adequately resolved and studied in previous studies as they used rudimentary lake models like a lumped (Piccolroaz et al., 2015) or a 1D lake model (Zhong et al., 2016). So, it is imperative to understand the influence of lake-atmosphere coupling on Great Lakes LST simulation and the mechanism behind it.
The two main objectives of this paper are to: (a) present a fully 3D two-way coupled lake-land-atmosphere modeling system for the Great Lakes using an RCM and a 3D lake model, and (b) assess the role of lake-atmosphere coupling on simulated summer LST. In this study, the Weather Research and Forecasting (WRF) model and an FVCOM-based 3D hydrodynamic lake and ice model are two-way coupled to each other. This fully 3D two-way coupled regional modeling system supplements the scarce modeling literature of RCMs with two-way coupled 3D lake models in the Great Lakes. Additionally, by performing a set of twin experiments where each experiment has a different coupling framework between WRF and FVCOM, this study aims to resolve and assess the impact of lake-atmosphere coupling on Lake Superior's summer LST in the two-way coupled model system. Lake Superior was chosen for this analysis because it has experienced an accelerated summer warming in recent decades (Austin & Colman, 2007;Zhong et al., 2016) and because Lake Superior, by virtue of being the deepest and largest among the Great Lakes, requires a comprehensive 3D lake model to accurately capture its hydrodynamics.
The remaining parts of this paper are organized as follows. In Sections 2 and 3, we describe the models and the model validation data used in this study along with the design of the numerical simulations. In Section 4, we validate the modeling results by comparing them to in-situ and satellite-based observation datasets. In Section 5, the role of lake-atmosphere coupling on Lake Superior's summer LST in the two-way coupled model system is presented, followed by discussion in Section 6 and conclusion and summary in Section 7.

Regional Climate Model
The RCM used in this study is the WRF model version 4.2.2 with the Advanced Research WRF (ARW) dynamic core (Skamarock & Klemp, 2008). The WRF physics opted for this study includes the Thompson microphysics scheme (Thompson et al., 2004(Thompson et al., , 2008, the Rapid Radiative Transfer Model for GCMs longwave and shortwave schemes (Iacono et al., 2008), the Yonsei University planetary boundary layer scheme (Hong & Lim, 2006;Noh et al., 2003) and the revised scheme of the surface layer formulation based on the fifth-generation Pennsylvania State University-National Center for Atmospheric Research Mesoscale Model (MM5) parameterization (Jiménez et al., 2012). Additionally, the land surface processes which include soil-vegetation physics are modeled by the Unified Noah land surface model within WRF (Chen & Dudhia, 2001).
The WRF domain for this study (Figure 1), centered at 85.0°W, 45.5°N, covers the entire Great Lakes region with 544 × 485 horizontal grid points at 4 km grid spacing. Vertically, the atmosphere is modeled with 50 stretched vertical levels topped at 50 hPa. The initial and boundary conditions for WRF to dynamically downscale are from the 3-hourly 0.25° fifth generation atmospheric reanalysis data (ERA5) (Hersbach et al., 2020) of the European Centre for Medium-Range Weather Forecasts (ECMWF). Here, no boundary nudging is applied so WRF is allowed to develop its own variability (e.g., spatial and internal variability) across the domain.
For the standalone WRF simulation used in this study (more details in Section 3), we used the satellite-derived LST and ice cover from the National Oceanic and Atmospheric Administration (NOAA) Great Lakes Environmental 4 of 24 Research Laboratory (GLERL) Great Lakes Surface Environmental Analysis (GLSEA) as the overlake surface boundary condition. GLSEA LST was found to be able to capture the very fine spatial features (1.3 km) of LST better than ERA5's SST data, leading to a better model performance in air temperature, and evaporation. See Wang et al. (2022) for more details about the standalone WRF model setup and difference in model performance when using different LST datasets as overlake surface boundary conditions.

Hydrodynamic Lake and Ice Model
In our two-way coupled modeling system, WRF v4.2.2 is interactively linked to an FVCOM-based 3D lake model to simulate the hydrodynamics, thermal dynamics, and ice dynamics of the Great Lakes. FVCOM is a prognostic, free-surface, 3D primitive equation coastal ocean circulation model that is numerically solved over an unstructured triangular grid using the finite-volume method (Chen et al., 2006(Chen et al., , 2013. As mentioned in Section 1, the Great Lakes have been successfully represented in 3D within regional climate modeling systems using FVCOM (e.g., Sun et al., 2020;Xue et al., 2017Xue et al., , 2022. Similar to Xue et al. (2022), the FVCOM variant used in this study is based on FVCOM version 4.1 without any nudging or other similar non-physical constraints so the lake hydrodynamic conditions freely interact with atmospheric conditions over the simulation period. WRF and FVCOM are run simultaneously with a two-way information exchange between them at 1-hr intervals using the OASIS3-MCT coupler (Craig et al., 2017). Here, the LST and ice cover are dynamically calculated by FVCOM and are provided to WRF as overlake surface boundary conditions. In turn, the atmospheric forcings required by FVCOM are dynamically calculated and provided by WRF.
The horizontal resolution of the FVCOM unstructured triangular grid ( Figure 2) ranges from ∼1 to 2 km near the coast to ∼2-4 km in the lake's offshore region. Vertically, the lakes are represented by 40 sigma layers to provide a vertical resolution ranging from <1 m in nearshore waters to ∼2-5 m in the lake's offshore regions. The Mellor-

of 24
Yamada level-2.5 (MY25) turbulence closure model (Mellor & Yamada, 1982) is used for simulating vertical mixing processes, including eddy viscosity and vertical diffusivities. The horizontal diffusivity, on the other hand, is calculated using the Smagorinsky numerical formulation (Smagorinsky, 1963). Finally, ice dynamics, ice thermal dynamics, and ice-water interaction processes are simulated by an unstructured-grid version of the Los Alamos Community Ice Code (CICE) embedded within the FVCOM framework (see Anderson et al. (2018) for more details).

Observation Data for Model Validation
The spatial pattern of seasonal air temperature and precipitation over the Great Lakes region from the standalone WRF and the WRF two-way coupled with FVCOM are evaluated against the Precipitation-Elevation Regressions on Independent Slopes Model (PRISM) data set developed by Daly et al. (1994Daly et al. ( , 1997Daly et al. ( , 2008. PRISM provides gridded temperature and precipitation data over the contiguous US at a spatial resolution of 1/8° latitude × 1/8° longitude. PRISM values are corrected for systematic elevation effects by using climate-elevation regression of weighted station data, with weights being assigned based on the station's physiography. Such correction is critical as elevation plays a major role in temperature and precipitation. Additionally, observation stations over high elevated regions (e.g., Appalachian, partially covered by our study domain) are preferentially located at lower elevations which leads to an underestimation of the true reading of precipitation, and overestimation of temperature, making corrections of elevation effects crucial. While PRISM has much higher spatial resolution in comparison to other widely used datasets such as the Global Air Temperature and Precipitation compiled by the University of Delaware (Willmott et al., 1998) and the Climatic Research Unit data (CRU) (Harris et al., 2020), it only covers the contiguous US. We, therefore, also use CRU at a grid spacing of 50 km (which covers the land of the entire globe) as the second source to evaluate the model performance of our standalone WRF and WRF two-way coupled with FVCOM. None of these observation datasets, however, have data over the Great Lakes. Thus, for air temperature, we also use the data from nine NOAA National Data Buoy Center (NDBC) scattered across the Great Lakes as shown in Figure 1.
The simulated latent and sensible heat fluxes over the lakes, which play a significant role in determining the LST of the lakes, are validated against the flux estimates derived from the NDBC buoy data. The fluxes were calculated using the Monin-Obukhov theory as described in Laird and Kristovich (2002) with slight modification to the Charnock coefficient to align it with the value used in WRF (see Wang et al. (2022) for more details).
Simulated LST are compared against GLSEA LST. Using satellite-derived imageries from the NOAA Advanced Very High Resolution Radar (AVHRR) and the Visible Infrared Imaging Radiometer Suite onboard the Suomi National Polar-Orbiting Partnership spacecraft (VIIRS S-NPP) and the NOAA-20 spacecraft, GLSEA produces one of the best available datasets to evaluate the spatial and temporal variability of the Great Lakes LST.

Design of Numerical Simulations
To achieve this study's objectives, two distinct model configurations ( Figure 3) are used to couple WRF and FVCOM. In the first configuration, WRF is run as a standalone model with the overlake surface boundary conditions (i.e., the LST and ice cover) prescribed from the daily gridded GLSEA data, as in Wang et al. (2022). GLSEA data is considered to be one of the best available representation of overlake surface boundary conditions for WRF to achieve an ideal WRF simulation; yet it is only available for historical simulations. FVCOM is then driven by the atmospheric outputs from the standalone WRF described above. This configuration is hereafter referred to as Observation-driven WRF-FVCOM (one-way) Coupling or OWFC in short. In the second configuration, WRF and the FVCOM model are run simultaneously with a two-way information exchange between them at 1-hr intervals using the OASIS3-MCT coupler. The LST and ice cover are dynamically calculated by FVCOM and are provided to WRF as overlake surface boundary conditions. Meanwhile, the atmospheric forcings required by FVCOM are dynamically calculated and provided by WRF. This configuration is hereafter referred to as WRF-FVCOM Two-way Coupling or WF2C in short. WF2C is designed to allow for as many variables as needed to be exchanged between WRF and FVCOM. A list of commonly used variables in different coupling scenarios in WF2C includes surface air temperature, surface air pressure, relative and specific humidities, total cloud cover, surface winds, downward shortwave radiation, downward longwave radiation, precipitation, evaporation, sensible and latent heat fluxes, LST, ice cover, and lake surface currents. The variables exchanged in this study under both OWFC and WF2C are illustrated in Figure 3. The difference between WF2C and OWFC is that the transfer of data between WRF and FVCOM is unidirectional in OWFC that is, FVCOM receives atmospheric forcing produced from standalone WRF simulations, and the overlake surface boundary conditions for WRF are prescribed using GLSEA data, rather than being dynamically calculated by FVCOM. On the contrary, WF2C eliminates the requirement of the observation data for the overlake boundary conditions in WRF. The LST and ice cover calculated by FVCOM, and the atmospheric variables calculated by WRF are allowed to freely interact with each other and evolve over time within WF2C. This advantage is pivotal as it sets the foundation for WF2C to be further developed to provide reliable future climate projections.
Since OWFC uses GLSEA data, which is one of the most accurate representations of LST and ice cover, to specify WRF boundary conditions while WF2C eliminates the requirement of the GLSEA data, validating WF2C's performance to be as comparable to that of OWFC would thereby attest to WF2C's ability to capture the regional climate and the lake-atmosphere dynamics of the Great Lakes. More importantly, the role of lake-atmosphere coupling on LST can be isolated and inferred from the discrepancy between the FVCOM simulated LSTs from WF2C and OWFC. Both WF2C and OWFC are initialized with the same initial condition for a continuous simulation period of May 2018 to December 2019, with May 2018 as the spinup period. Aimed for understanding the impact of lake-atmosphere coupling on summer LST, when the lake thermal structure is complex with stratification, the validation and analysis for this study are focused on the summer (June, July, August (JJA)) and fall (September, October, November (SON)) seasons of 2018 and 2019.

Air Temperature and Precipitation
As mentioned in Section 3, as OWFC relies on GLSEA data for an accurate representation of LST and ice cover to drive WRF, whereas WF2C does not, we expect that if the two-way coupling between WRF and FVCOM performs well, then OWFC and WF2C will have similar performance in reproducing the air temperature and precipitation over the Great Lakes region when compared to CRU, PRISM and NDBC buoys (Figures 4, 6, and 7). The comparison of spatial distribution of seasonal air temperature from PRISM, CRU, OWFC, and WF2C are shown in Figure 4. CRU, which has a coarser spatial resolution than PRISM, produces a similar spatial distribution of air temperature to PRISM for the US Great Lakes region. As shown in Figure 4, overland, the spatial correlations of WF2C with CRU and PRISM, which are 0.943 and 0.951 for JJA, respectively, and 0.975 and 0.980 for SON, respectively, are similar or slightly higher than the correlation between OWFC and the observation datasets. WF2C captures the meridional gradient in air temperature well with higher air temperature simulated for lower latitudes. Higher air temperature observed in the southwestern portion during JJA is also replicated by WF2C. The root mean square difference (RMSD) of WF2C with CRU and PRISM are relatively small with 1.565°C and 0.832°C for JJA respectively and 1.527°C and 0.929°C for SON respectively. These RMSDs of WF2C are similar or slightly smaller than the RMSDs of OWFC. Also, in general, the correlations are higher and RMSDs are smaller when comparing the models with PRISM than with CRU, indicating the importance of high spatial resolution observation datasets when evaluating high resolution models.
Over the lakes, both WF2C and OWFC have very high spatial correlation to each other but WF2C simulates noticeably warmer air than OWFC during JJA, especially over Lake Superior ( Figure 5). The air over the central basin of Lake Superior is warmer by up to 3°C in WF2C relative to OWFC. The warmer overlake air also affects Lake Superior's summer LST which is discussed in Section 4.3 and Section 5. In comparison to the overlake air temperature from the buoys across the Great Lakes, WF2C and OWFC both reproduce the air temperature remarkably well ( Figure 6). They perform similarly to each other with both WF2C and OWFC having a correlation >0.945 and an RMSD <1.5°C with buoy observations, except for the month of July in Lake Superior (Figures 6a, 6d, and 6f), where WF2C overestimates the air temperature.
Although the primary focus of this study is on summertime, and both WF2C and OWFC have been successfully calibrated and validated for summer and fall, it's worth noting that the current WRF configuration exhibits a relatively larger cold bias for near-surface air temperature during winter and spring, reaching up to 5°C over some locations. Employing the Noah-MP land surface model (Niu et al., 2011) can partially mitigate this cold bias in air temperature during these seasons, but it tends to produce a larger warm bias in other seasons. As such, future investigations will be necessary to address the limitations of the WRF model when simulating colder seasons. Figure 7 compares the summer and fall precipitation patterns over the Great Lakes region between the observation datasets and the models. Unlike in air temperature, PRISM and CRU differ noticeably with PRISM unsurprisingly having finer spatial variations due to its higher spatial resolution. Both OWFC and WF2C successfully reproduce PRISM's finer spatial variations in precipitation that are absent in CRU such as the higher precipitation upwind of the Great Lakes (e.g., in Wisconsin) likely caused by mesoscale convective systems and downwind of Great Lakes (e.g., southeast of Lake Erie) likely caused by isolated deep convections (Wang et al., 2022). As such, this serves as an excellent illustration of why it is preferable to use high-resolution observation datasets like PRISM for validating the models. Upon closer examination, we see the WF2C performs slightly better than OWFC when validating against PRISM. Compared to OWFC, WF2C exhibits a marginally smaller RMSD when compared with CRU and PRISM with 1.077 mm/day and 1.324 mm/day for JJA respectively and 1.024 mm/day and 1.217 mm/day for SON respectively. The spatial correlations of WF2C with CRU and PRISM are 0.461 and 0.362 for JJA respectively and 0.523 and 0.568 for SON respectively, which are generally slightly higher than the correlation of OWFC with CRU and PRISM (Figure 7).

Latent and Sensible Heat Fluxes
Figure 8 presents a temporal and statistical comparison of latent and sensible heat fluxes between the buoys, WF2C, and OWFC. Both WF2C and OWFC effectively reproduce the seasonality of the fluxes. They capture not only the overall magnitude of the fluctuations, but also the rapid changes in sensible heat on daily to weekly scales in the fall (Figure 8a-8j). Compared to OWFC, statistically, WF2C has a noticeably higher correlation with observations as well as lower RMSD (Figures 8k and 8l). In particular, for Buoy 45001 in Lake Superior, during July, OWFC overestimates a large outgoing latent and sensible heat fluxes while the buoy data suggests that the fluxes are incoming fluxes of relatively smaller magnitudes. WF2C, in contrast, aligns more closely with the buoy data. For latent heat flux, the average correlation between observation and WF2C is higher than between observation and OWFC by 0.05 while the average RMSD of WF2C is smaller than that of OWFC by 2.85 W/m 2 . For sensible heat flux, the average correlation between observation and WF2C is higher than that between observation and OWFC by 0.06 while the average RMSD of WF2C is smaller than the average RMSD of OWFC by 2.06 W/m 2 .

Lake Surface Temperature
The spatiotemporal comparison of LST between GLSEA, OWFC, and WF2C is shown in Figure 9. During June to November, the Great Lakes undergo an increase and a subsequent decrease in LST with distinct spatial differences in the meridional direction. As the northernmost and the deepest lake, Lake Superior is always the coolest lake while Lake Erie, being the southernmost and shallowest lake, is the warmest lake. This warming gradient is particularly noticeable in July as Superior maintains a lakewide average LST of ∼13°C while Lake Erie uniformly warms to a lakewide average of ∼24°C. Both WF2C and OWFC capture this inter-lake variation very well, although both produce a noticeable warm bias for Lake Superior in July. The models are also able to capture the spatial heterogeneity within each lake such as the distinct north-south gradient in Lakes Michigan during July-August (Figures 9e-9k) and the distinct east-west gradient in Lake Erie during September (Figures 9m-9o).
In addition to resolving the spatial variability of monthly LST, WF2C and OWFC also track the magnitude and the temporal evolution of the lakewide average LST well (Figure 10). Both WF2C and OWFC have good  correlation (>0.960) and RMSD (<2°C) when compared with GLSEA. Looking at the previous studies that used 1D (Bennington et al., 2014;Notaro et al., 2015) and 3D lake models (Bai et al., 2013;Xue et al., 2017Xue et al., , 2022, it is clear that such a close tracking of LST, particularly the spring-early summer warming and the summer peaks, is only achievable through the use of 3D lake models. The models, however, do have some biases including a warm and cold bias during July-September of Lake Superior (∼2°C) and a persistent cold bias in Lake Ontario (∼2-3°C). Such LST biases can be expected as modeling the physical processes of deep lakes is challenging (Bai et al., 2013;Bennington et al., 2014;Xiao et al., 2016;Xue et al., 2017) and Lake Superior and Lake Ontario are the deepest and second deepest lakes among the Great Lakes in terms of average depth, respectively. Nevertheless, when compared to a more constrained, atmospheric data-driven driven FVCOM simulation such as the standalone FVCOM driven by reanalysis data in Bai et al. (2013), both WF2C and OWFC perform admirably. For example, in Lake Superior and Erie, both WF2C and OWFC produce a lower lakewide average LST RMSD than Bai   average LST RMSD of 1.14°C while both WFC and OWFC have just slightly higher RMSDs of 1.80°C and 1.91°C respectively.
The performance of WF2C and OWFC in reproducing LST is quite similar, with WF2C generally achieving a similar, if not slightly lower, RMSD than OWFC for the lakewide average LST (by up to 0.117°C) for all lakes, with the exception of Lake Superior. The apparent discrepancies in Lake Superior between WF2C and OWFC during July and August is noteworthy (Figures 9h, 9l, and 10a). For July, while both OWFC and WF2C produce a warm bias for the lakewide average LST, WF2C has a higher bias of ∼2°C and OWFC has a lower bias of ∼1°C. It is interesting to note that even though OWFC overestimates a large outgoing latent and sensible flux (Figures 8a  and 8b)-which would favor a cooling of LST-it still produces a warm LST bias in July. This suggests that changes in other variables within OWFC, such as radiation fluxes and lake mixing, are compensating for the outgoing latent and sensible fluxes, preventing the LST from dipping below the observed value. Similarly, even though WF2C's latent and sensible heat fluxes closely track the observation, it still has a ∼2°C bias, implying that other variables must be playing a role in warming the simulated LST. For August, WF2C accurately captures the peak LST magnitude while OWFC underestimates the peak by ∼2°C. This discrepancy is likely due to their respective simulated antecedent (July) LST, as well as the synchronic surface heat fluxes and lake hydrodynamic conditions in August. The difference between WF2C and OWFC fundamentally stems from the presence and absence of information exchange between FVCOM and WRF that is, lake-atmosphere coupling. It is therefore possible to infer and examine the impact that lake-atmosphere coupling has on Lake Superior's summer LST through our twin experiment of WF2C and OWFC (discussed in Section 5).

Impact of Coupled Lake-Atmosphere Dynamics in Summer LST Simulation
Understanding the impact of coupled lake-atmosphere dynamics on the simulation of the Great Lakes LST is critical not only for explaining observed historical warming but also for ensuring the accuracy of future LST projections, which inevitably rely on a coupled lake-atmosphere system. Current understanding is hampered by the prevalent use of 1D lake models in regional climate modeling, which results in an insufficient representation of coupled lake-atmosphere dynamics and leads to a drift in simulated LST from its true state. The two-way coupled WRF-FVCOM system offers us the opportunity to examine how lake-atmosphere coupling could impact the Great Lakes LST by deriving insights from the discrepancies between WF2C and OWFC simulations. In this section, we focus on the role of lake-atmosphere coupling in influencing the summer LST simulation of Lake Superior, the largest and deepest one among the five lakes, which contains more than 50% of total water mass of the Great Lakes. Lake Superior's hydrodynamic summer is typically considered to be from July-September, a period during which the lake exhibits the warmest surface temperature. Hence, our analysis in this section focuses on the period from late spring to the end of the hydrodynamic summer, that is, from June to September.
Our experiment shows that coupled lake-atmosphere dynamics increases the simulated summer LST (i.e., the differences between WF2C and OWFC) as shown in Figure 11. LST from WF2C is consistently higher than LST from OWFC during June-September. The trend and magnitude of the LST difference between WF2C and OWFC are very similar in both years, suggesting it is a consistent pattern due to lake-atmosphere coupling rather than being an isolated episodic event. The increase in LST commences from mid-June, reaching a peak in late July or early August with a magnitude of approximately +1.3°C. The increase in LST is then more or less sustained until mid-August after which it suddenly drops and gradually levels out at around +0.3°C.
The lake-atmosphere coupling elevates the summer LST by modifying the surface heat fluxes going into or out of the lake. The surface heat fluxes responsible for LST changes are upward longwave radiation (ULW), downward longwave radiation (DLW), sensible heat flux (SH), latent heat flux (LH) and net shortwave radiation (NSW). The Net Heat (aggregate of ULW, DLW, SH, LH and NSW) that the lake receives is the primary driving factor for lake warming. The magnitude of the surface heat fluxes in both WF2C and OWFC for June 2018 are compared in Figure 12. By examining the differences between WF2C and OWFC for June 2018, we see that the lake-atmosphere coupling results in a higher Net Heat into the lake (by 11.02 W/m 2 ). This increase in Net Heat is primarily due to the increase in the SH entering the lake (by 5.01 W/m 2 ) and decreases in the LH exiting the lake (by 6.41 W/m 2 ). Figure 13a summarizes these differences in surface heat fluxes between WF2C and OWFC for June 2018. Figures 13b-13h similarly encapsulate the differences for other months. Figure 13 shows that, due to the lake-atmosphere coupling, the lake gains more energy (positive Net Heat change) in June and July and gains less energy (negative Net Heat change) in August and September. The most significant heat gain occurs in July, with 22.07 W/m 2 in 2018 and 23.33 W/m 2 in 2019. The largest loss is in September, with 16.26 W/m 2 in 2018 and 29.21 W/m 2 in 2019. The energy gain (loss) is primarily due to the decrease (increase) of outgoing latent heat and increase (decrease) in the incoming sensible heat. Importantly, note that the changes in Net Heat between WF2C and OWFC and the changes in LST between WF2C and OWFC are not synchronous due to the lake seasonal mixing, as depicted in Figure 14.
In the unstratified waters of June, even though more heat (positive Net Heat change) is directed into the lake in WF2C, the surplus heat disperses throughout the water column, leading to a very minor increase in LST when compared to OWFC. Then in July, as shown in Figure 11, the LST in WF2C undergoes a rapid increase relative to OWFC. This rapid warming is due to two reasons. First, the surplus heat is the largest in July and second, the water is stratified in July, which creates an ideal condition for the warming to be restricted just to the surface layer of the water. Later in August, less heat (negative Net Heat change) enters the lake in WF2C compared to OWFC, which curtails the rapid LST increase seen in July. However, the reduced heat input in WF2C, combined with the fully stratified water, is sufficient to sustain (but not further increase) the LST difference between WF2C and OWFC established in July. Finally, in September, the combination of weakening stratification and further reduction of the incoming heat in WF2C relative to OWFC leads to a diminishing LST difference between the two simulations, culminating in a convergence of the LSTs from the two simulations ( Figure 11). It is intriguing to observe the slightly warmer deeper waters of OWFC during August and September shown in Figure 14. This occurs, in part because, in the WF2C model, a higher LST along with stronger stratification slows down the vertical water mixing, thereby keeping the deeper water cooler during the July-August period, while keeping the upper layer warmer. Conversely, in the OWFC model, the relatively weaker stratification facilitates an easier distribution of heat throughout the water column. This results in the deeper waters of OWFC being slightly warmer than those in WF2C, with a correspondingly cooler upper layer.  Table S1 in Supporting Information S1 for a tabular format. Figure 13. Differences between WF2C and OWFC in surface heat fluxes from June-September 2018 (left panels) and June-September 2019 (right panels). The arrows pointing down (yellow) are surface heat fluxes that put heat into the lake. The arrows pointing up (red) are surface heat fluxes that take heat from the lake. Negative values mean the flux has a lower magnitude in WF2C than OWFC. Positive values mean the flux has a higher magnitude in WF2C than OWFC. See Tables S1 and S2 in Supporting Information S1 for a tabular format.
As discussed above, lake-atmosphere coupling affects the modeled LST by modifying surface heat fluxes between WF2C and OWFC. This is realized through lake-atmosphere coupling influencing various meteorological state variables interacting with LST. The impact of lake-atmosphere coupling on surface heat fluxes and the associated meteorological state variables are demonstrated in Figure 15 for 2018 and Figure 16 for 2019. For instance, the changes in SH are mainly due to the changes in the difference between air temperature and LST. In June and July, the air temperature increases more rapidly than the LST. This leads to a larger increase in the difference between air temperature and LST, resulting in an increase in SH in the 2 months when comparing WF2C and OWFC simulations (Figures 15j, 15l, 16j, and 16l). Similarly, the changes in LH (Figures 15p and 16p) are due to the changes in the difference between saturated and actual specific humidity (Figures 15n and 16n). Although wind speed is also one of the factors that affects sensible and latent heat flux, the changes in wind speed averaged over the entire lake are relatively small (Figures 15v and 16v).
These causal relationships of sensible and latent heat flux with the meteorological state variables are represented by the bulk transfer equation of heat fluxes (Arya, 2001). The ULW and DLW are also largely controlled by the state of LST and air temperature, as expressed by the Stefan-Boltzmann law of radiation (Figures 15b, 15f, 16b, and 16f). These interactive processes further impact atmospheric conditions such as cloud formation and NSW. However, it is important to note that more comprehensive analysis is needed to disentangle the detailed processes that led to these changes. This includes examining both large-scale and local forcings, along with the interactions among various atmospheric state variables.

Sensitivity of Model LST to Varying Coupling Approaches
In the previous section, we showed that the warm bias in Lake Superior in July is larger in WF2C than in OWFC. While there are a number of factors in the model configurations that affects model performance-such as grid resolution, WRF domain and lateral boundary forcing, WRF parameterization schemes (for microphysics, cumulus convection, planetary boundary layer, radiation, and land surface), and FVCOM turbulent schemes-this is not the main focus of this study. This study concentrates on the impact of lake-atmosphere coupling on LST, as reflected in the difference between WF2C and OWFC simulations under identical model configuration except with and without two-way coupling. With that in mind, we wanted to ensure that such a difference was not caused by any potential inconsistencies in the calculation of fluxes in WRF and FVCOM. Specifically, in the state-variable-based coupling described above, WRF and FVCOM directly exchange state variables and calculate fluxes within each model separately; meteorological variables such as wind, air temperature, cloud coverage, and relative humidity are transferred from WRF to FVCOM. This method, widely used in coupled ocean-atmosphere or lake-atmosphere coupling, can provide more accurate results by leveraging the different grid resolutions in atmospheric and hydrodynamic models (Sun et al., 2020;Warner et al., 2010;Xue et al., 2015Xue et al., , 2017, given both models have compatible formulations for computing the fluxes. Within FVCOM, the Coupled Ocean-Atmosphere  First column: magnitude of lake surface temperature (LST) (a), near surface air temperature (T2) over the lake (e), near surface air temperature over the lake minus LST (i), saturated specific humidity minus near-surface specific humidity over the lake (m), cloud cover over the lake (q) and wind speed over the lake (u) for 2018. Second column: difference between WF2C and OWFC for the adjacent first column panel's atmospheric variable. Third column: magnitude of ULW (c), DLW (g), SH (k), LH (o), NSW (s) and Net Heat (w) for 2018. Positive fluxes represent surface heat fluxes that put heat into the lake while negative fluxes represent surface heat fluxes that take heat from the lake. Fourth column: difference between WF2C and OWFC for the adjacent third column panel's surface heat fluxes. See Table S1 in Supporting Information S1 for exact numbers. Figure 16. First column: magnitude of lake surface temperature (LST) (a), near surface air temperature (T2) over the lake (e), near surface air temperature over the lake minus LST (i), saturated specific humidity minus near-surface specific humidity over the lake (m), cloud cover over the lake (q) and wind speed over the lake (u) for 2019. Second column: difference between WF2C and OWFC for the adjacent first column panel's atmospheric variable. Third column: magnitude of ULW (c), DLW (g), SH (k), LH (o), NSW (s) and Net Heat (w) for 2019. Positive fluxes represent surface heat fluxes that put heat into the lake while negative fluxes represent surface heat fluxes that take heat from the lake. Fourth column: difference between WF2C and OWFC for the adjacent third column panel's surface heat fluxes. See Table S2 in Supporting Information S1 for exact numbers.
Response Experiment (COARE3) (Charusombat et al., 2018;Edson et al., 2014;Fairall et al., 2003) scheme is used to calculate surface fluxes based on these atmospheric variables and the LST. Conversely, lake temperature and ice cover from FVCOM serve as WRF's over-lake boundary conditions, and WRF uses the revised MM5 surface layer scheme to calculate surface fluxes (Jiménez et al., 2012). While both schemes are based on the Monin-Obukhov similarity theory, and the revised MM5 surface layer scheme indeed adopts the COARE3 similarity function for unstable atmospheric conditions, it is important to ensure that both schemes compute compatible fluxes for consistency. In other words, we need to ensure that the warmer LST in WF2C in Lake Superior in July is not due to the use of different surface flux schemes in FVCOM and WRF. Therefore, we implemented an alternative flux-based coupling approach for cross-validation Xue et al., 2014Xue et al., , 2020. In the flux-based coupling, WRF and FVCOM directly exchange surface fluxes of energy, momentum, and mass across the lake-atmosphere interface. Specifically, FVCOM provides WRF with LST and ice cover as its overlake boundary conditions. Meanwhile, FVCOM directly receives calculated fluxes from WRF as driving forces, instead of calculating surface fluxes within FVCOM. These fluxes include heat/radiation fluxes (sensible, latent, longwave, and shortwave radiation) and momentum flux (wind stress). This flux-based coupling approach ensures strict consistency in the fluxes exchanged between the lake and atmosphere within both the atmospheric model and lakes, as they are all calculated within WRF. However, it does not allow WRF and FVCOM to leverage their varying spatial resolutions to calculate surface fluxes.
In our case, the two coupling approaches-flux-based and state-variable-based-exhibit highly consistent results when the grid resolutions of WRF and FVCOM are similar (WRF: 4 km and FVCOM: 1-4 km). This comparison addressed our concern about any potential inconsistency arising from using the revised MM5 surface layer scheme in WRF and the COARE3 scheme in FVCOM. As shown in Figure 17, there is only a 2% difference (∼0.16°C) in the simulated LSTs between the flux-based coupling and state-variable-based coupling in our case, with the warming patterns being nearly identical. Therefore, the differences in warming between WF2C and OWFC indeed stems from the two-way coupling of WRF and FVCOM and is not related to the different surface flux schemes in FVCOM and WRF.

Computational Cost
The computational effort differs between running WF2C and OWFC, with WF2C requiring approximately 1.1-1.3 times more computational time than OWFC. This is due to the end-to-end coupling process in WF2C. The increase in computational time stems from two major factors: (a) the extra time required for data exchange at each coupling time step, and (b) the inconsistent computing time required by the two coupled components. That is, the slower model determines the overall performance, as both models must reach the same simulation step before data exchange can occur, causing the faster model to idle while waiting for the slower one. This difference in computation time could be significant, especially for multi-decadal simulations. Given the comparable Figure 17. Comparison of simulated lake surface temperature using flux-based coupling and state-variable-based coupling for Lake Superior.
performance of WF2C and OWFC, OWFC could still be a viable option when one is only interested in historical simulations. However, OWFC is not suitable for future projections due to its inherent reliance on a robust observational LST data set, which is required by the atmospheric model (WRF) as the lower boundary condition.

WF2C's Potential for Long-Term Climate Simulation
We acknowledge that more extended simulations are crucial for a comprehensive evaluation of the WF2C's performance in long-term climate simulation. However, it's important to highlight that our 2-year simulation yielded valuable insights, particularly in the context of long-term climate modeling. Our 2-year simulation was designed to investigate whether the observed warming bias in the WF2C would progressively intensify, thereby causing the model to gradually deviate from real-world conditions. If this scenario were to manifest, WF2C would fail to qualify for further application in long-term, climate-scale simulations. In contrast, we observed that the warm bias emerged annually, following a consistent pattern, and crucially, this bias subsided and lessened during the autumn months in the lake-atmosphere coupled simulation. This implies that the warm bias does not perpetually accumulate over years; thereby, providing us with confidence in WF2C's potential for further development for long-term simulations and future projections.

Conclusion and Summary
In this study, we developed a fully coupled modeling system (WF2C) with WRF and FVCOM to accurately capture the Great Lakes and its two-way interaction with the atmosphere. We also developed an observation-driven WRF-FVCOM (one-way) Coupling (OWFC) in which we first ran a standalone WRF that uses satellite-derived LST and ice cover from GLSEA as overlake surface boundary conditions and then used the standalone WRF output to drive FVCOM to simulate the 3D status of the lakes. The fundamental distinction between WF2C and OWFC simulation lies in the fact that in WF2C, the overlake surface boundary condition for WRF is dynamically calculated by FVCOM, while in OWFC, the overlake surface boundary condition for WRF needs to be prescribed using observational data of LST and ice cover (GLSEA). WF2C and OWFC were validated against in-situ and satellite-based observation datasets, with both showing strong performance in reproducing the historical conditions of the Great Lakes during summer and fall. OWFC simulations, due to their dependence on satellite-derived LST and ice cover for defining overlake surface boundary conditions, are suited only for historical analysis. On the other hand, W2FC can potentially be utilized for future climate projections as it allows the lake hydrodynamic conditions and atmospheric conditions to evolve internally, with the atmosphere and lake and freely interacting with each other over the entire course of climate simulation. The fact that WF2C has achieved similar performance to the OWFC while eliminating the requirement of the GLSEA data is exciting because it sets the foundation for future studies to explore WF2C's potential to produce reliable climate projections and accurately simulate the lake-atmosphere dynamics of the Great Lakes under future climate conditions that could differ significantly from climate condition represented in the observational record.
We also used the WF2C and OWFC simulations to analyze the role of lake-atmosphere two-way coupling in affecting the simulated components of the summertime lake-atmosphere surface heat fluxes over Lake Superior. This analysis is vital not only for accurately simulating LST, but also for effectively evaluating, understanding, and documenting the behavior of the newly developed model. Acknowledging both its strengths and limitations is crucial to the model's continual development. Our study highlights that in any coupled-model development, the two-way coupling process inherently permits free interactions between different model components. This allows for unimpeded propagation of dynamic processes and sometimes also leads to amplification of perturbations and errors within the coupled system. The warm biases in this case, are an example of how errors can propagate and amplify, but they could also manifest as cold biases or other formats, depending on model configurations. It is essential that these phenomena are not only accurately documented, but also thoroughly understood by the model development community.
By inferring the impact of coupled lake-atmosphere dynamics from the difference between WF2C and OWFC, our results show that lake-atmosphere coupling in WF2C can lead to higher LST during the summer period through modifying surface heat fluxes and influencing various meteorological state variables interacting with LSTs. Notably, although OWFC seems to have a smaller warm bias in simulating July LST than WF2C does, 10.1029/2023MS003620 22 of 24 OWFC actually overestimates the outgoing latent and sensible heat fluxes when compared to observation. Furthermore, WF2C has a larger warm bias of July LST and overlake air temperature than OWFC despite having a better reproduction of latent and sensible heat fluxes than OWFC, as discussed in Section 4. The overestimation of LST can be caused by a combination of various factors including the atmospheric (Wang et al., 2022) and hydrodynamic processes (Ye et al., 2019) simulated in each model as well as the lake-atmosphere coupling.
Comparisons in the simulated surface heat fluxes in WF2C and OWFC clearly reveal new processes are activated in the WF2C simulation, underscoring the importance of coupling for accurate prediction of lake-atmosphere heat exchange. These mechanisms are currently not well understood and warrant more focused investigation in the future Great Lakes studies. These results highlight the complexity of the coupled lake-atmosphere system and sheds light on the importance of process-oriented studies to reveal key processes that influence the model state and allow us to improve regional climate simulations. This study identified lake-atmosphere coupling as a necessary feature of Great Lakes hydrodynamic simulations, as it provides reasonably accurate estimates of summer LST and surface heat fluxes even for simulation unconstrained by observations. Lastly, while our study focused on summertime lake-atmosphere interaction and their impacts on LST, it is important to note that during the colder months of the year, snow and ice are also a major part of the lake-atmosphere interaction. For example, ice cover over the lakes acts as a physical barrier between the lake and the atmosphere, thus affecting the surface heat fluxes (Fujisaki-Manome et al., 2020;Vavrus et al., 2013). These new variables, which come into play during the colder seasons, add another layer of complexity to the coupled system. This complexity is exemplified by our WRF configuration which produces a relatively larger cold bias during winter and spring while performing remarkably well during summer and fall. So future studies should carefully consider snow and ice cover to provide an annual picture of the importance of lake-atmosphere coupling on the Great Lakes hydrodynamic and regional climate simulations.