dataesr/scientific-paragraphs-categorization huggingface dataset:
Materials and methods This section first describes the study domain (Sect. 2. 1), datasets used to set up and force models (Sect. 2. 2), and the common routing product produced for this study and used by most of the models (Sect. 2. 3). A brief description of the 13 participating models can be found in Sect. 2. 4. The datasets used to calibrate and validate the models regarding streamflow and the datasets to evaluate the models beyond streamflow are introduced in Sects. 2. 5 and 2. 6, respectively, followed by Sect. 2. 7, which defines the metrics to determine the models' ability to simulate the variables provided in these datasets using two different approaches, depending on which spatial aggregation is applied (Sect. 2. 8). A multi-objective analysis described in Sect. 2. 9 is performed to augment the previous single objective analyses that evaluated the performance across the four variables independently.
5 %). These estimates are derived using the vertically averaged soil data based on the entire soil column up to 2. 3 m from the Global Soil Dataset for Earth System Models (GSDE;  Shangguan et al. , 2014)  dataset, and soil classes are based on the United States Department of Agriculture's (USDA) classification. The elevation of the Great lakes watershed (including the Ottawa River watershed) ranges from 17 to 1095 m. The mean elevation is 270 m, and the median is 247 m. These estimates are derived using the HydroSHEDS 90 m (3 arcsec) digital elevation model (DEM) data  (Lehner et al. , 2008). These datasets -NALCMS, GSDE, and Hy-droSHEDS -summarize the common datasets used by all partners in this study to setup models. More details about the datasets are provided in Sect. 2. 2.
Meteorologic forcings and geophysical datasets All models contributing to this study used the same set of meteorologic forcings and geophysical datasets to set up and run their models. Additional datasets required to set up specific models had to be reported to the team and made available to all collaborators in order to make sure everyone would have the chance to use that additional information as well. The common datasets used here were determined after the preceding Great Lakes Runoff Intercomparison project for Lake Erie (GRIP-E;  Mai et al. , 2021)  in which contributors were allowed to use datasets of their preference. The team for the Great Lakes project (Great Lakes Runoff Intercomparison Project Phase 4: the Great Lakes -GRIP-GL) evaluated all datasets used by all models in GRIP-E and decided together which single dataset to commonly use in GRIP-GL.
This led to the situation in which most models had to be set up again as, for all models, at least one dataset differed from what was used in GRIP-E and what was decided to be used in GRIP-GL. The common datasets are briefly described in the following. The Regional Deterministic Reanalysis System v2 (RDRS-v2) was employed as meteorologic forcing  (Gasset et al. , 2021). RDRS-v2 is an hourly and 10 km by 10 km meteorologic forcing dataset covering North America (Fig. 2a ). Forcing variables such as precipitation, temperature, humidity, and wind speed, among others, are available from January 2000 until December 2017 through the Canadian Surface Prediction Archive (CaSPAr;  Mai et al. , 2020b). The dataset has been provided to all the participating groups in two stages. First, only data for the calibration period (January 2000 to December 2010) were shared.
Second, after finalizing the calibration of all models, the data were shared for the validation period (January 2011 until December 2017). This allowed for a true blind (temporal) validation of the models. This was possible since the RDRS-v2 dataset was not pub-licly available when the GRIP-GL project started, and the participating groups (except the project leads) did not have access to the data. The forcings were aggregated according to the needs of each model. For example, lumped models used aggregated basin-wise daily precipitation and minimum and maximum daily temperature instead of the native gridded, hourly inputs. The HydroSHEDS dataset  (Lehner et al. , 2008)  was used as the common digital elevation model (DEM) for the project. It has a 3 arcsec resolution which corresponds to about 90 m at the Equator.
Since the upscaled HydroSHEDS DEMs with 15 arcsec (500 m) and 30 arcsec (1 km) are consistent with the best-resolution dataset, the collaborators were allowed to pick the resolution most appropriate for their setup. The DEM was then postprocessed to the specific needs of each model. For example, some lumped models required the average elevation of each watershed as an input. The data were downloaded from the HydroSHEDS website (https://www. hydrosheds. org, last access: 6 July 2022), cropped to the study domain, and provided to the collaborators. The Global Soil Dataset for Earth System Models (GSDE;  Shangguan et al. , 2014)  was used as the common soil dataset for all models in this study. This dataset, with a 30 arcsec spatial resolution (approximately 1 km at the Equator), contains eight layers of soil up to a depth of 2.
3 m. The data were downloaded directly from the website mentioned in the  Shangguan et al. (2014)  publication and preprocessed for the needs of the collaborators and models. For example, the data were regridded to the grid of the meteorologic RDRS-v2 forcings, as this was the grid for which most distributed models were set up. Some other models required the aggregation of soil properties to specific (fixed) layers different to those distributed. The project leads also converted the soil textures provided to soil classes as this was required by a few models. All these data products were converted to a standardized NetCDF format and shared with the collaborators. As a common land cover dataset for all models, the North American Land Change Monitoring System (NAL-CMS) product including 19 land cover classes for North America was used.
The dataset has a 30 m by 30 m resolution and is based on Landsat imagery from 2010 for Mexico and Canada and from 2011 for the USA. The data can be downloaded from a website  (NACLMS, 2017). The data were downloaded, cropped to the study domain, and saved in a common TIFF file format. Furthermore, model-specific postprocessing of the data included regridding to model grids and aggregation of the land cover classes to the eight classes common in MODIS. Those datasets were saved in standard NetCDF file formats and provided to the collaborators.
0;  Han, 2021) , which is a set of Pythonbased geographic information system (GIS) tools for supporting vector-based watershed delineation, discretization, and simplification and incorporating lakes and reservoirs into the network. This toolbox has been successfully applied for producing the Pan-Canadian  (Han et al. , 2020)  and North American routing product  (Han et al. , 2021b). More detailed information, a user manual, and tutorial examples can be found on the BasinMaker website  (Han et al. , 2021a). In this study, we used the HydroSHEDS DEM, flow direction, and flow accumulation at the 3 arcsec resolution for BasinMaker to delineate the watersheds and discretize them into subbasins. The HydroLAKES database, which provides a digital map of lakes globally  (Messager et al. , 2016) , is used in the watershed delineation as well.
The Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1) Version 6 data (USGS, 2019) is used in BasinMaker to estimate the Manning's n coefficient for the floodplains. The global river channel bankfull width database developed by  Andreadis et al. (2013)  is used to calculate geometry characteristics of river cross sections. This routing product is compatible with other hydrological models (i. e. , runoff-generating models at any resolution, either gridded or vector based) when the Raven hydrologic modeling framework  (Craig et al. , 2020)  is used as the routing module. The routing product contains adequate information for hydrologic routing simulation, including the rout-ing network topology, attributes of watersheds, lakes, and rivers, and initial estimates for routing parameters (e.
, 2021). The figure shows the mean annual precipitation derived for the 18 years during which the dataset is available and used in this study  (2000) (2001) (2002) (2003) (2004) (2005) (2006) (2007) (2008) (2009) (2010) (2011) (2012) (2013) (2014) (2015) (2016) (2017). The forcing inputs are preprocessed as needed by the models, e. g. , aggregated to basin averages or aggregated to daily values. (b) Actual evapotranspiration estimates and (c) surface soil moisture from the Global Land Evaporation Amsterdam Model  (GLEAM v3. 5b; 25 km; daily; Martens et al. , 2017)  and (d) snow water equivalent estimates from the ERA5-Land dataset (Muñoz Sabater, 2019; 10 km; daily) have been used to evaluate the calibrated model setups. All datasets have been cropped to the study domain (black line). The actual evapotranspiration and surface soil moisture (panels b and d) are not available over the lake area (missing grid cells).
(2018). The model was set up using the NeuralHydrology Python library  (Kratzert et al. , 2022 ). An LSTM is a deep-learning model that learns to model relevant hydrologic processes purely from data. The LSTM setup used in this study is similar to that from  Kratzert et al. (2019c)  and has been successfully applied for streamflow prediction in a variety of studies (e. g. ,  Klotz et al. , 2021; Gauch et al. , 2021; Lees et al. , 2021). The model inputs are nine basin-averaged daily meteorologic forcing data (precipitation, minimum and maximum temperature, u/v components of wind, etc. ), nine static scalar attributes derived for each basin from these forcings for the calibration period (mean daily precipitation, fraction of dry days, etc. ), 10 attributes derived from the common land cover dataset (fraction of wetlands, cropland, etc.
), six attributes derived from the common soil database (sand, silt, and clay content, etc. ), and five attributes derived from the common digital elevation map (mean elevation, mean slope, etc. ). ). The Supplement contains the full list of inputs, and a more detailed overview of the training procedure (Sect. S3. 1). Unlike process-based hydrologic models, the internal states of an LSTM do not have a direct semantic interpretation but are learned by the model from scratch. In this study, the LSTM was only trained to predict streamflow, which is why it did not participate in the comparison of additional variables in this study. Note, however, that even though the model does not explicitly model additional physical states, it is nevertheless possible -to some degree -to extract such information from the internal states with the help of a small amount of additional data  (Lees et al.
g. , used by the U. S. Army Corps of Engineers in long-term water balance monitoring and seasonal water level forecasting). The model will be referred to as LBRM-CC-lumped hereafter. Inputs to LBRM-CC-lumped include daily areally averaged (lumped) precipitation, maximum temperature, and minimum temperature, as well as the contributing watershed area. All these inputs are based on the common dataset described in Sect. 2. 2. For each basin, LBRM-CC-lumped's nine parameters were calibrated against the Kling-Gupta efficiency (KGE;  Gupta et al. , 2009)  using a dynamically dimensioned search (DDS) algorithm  (Tolson and Shoemaker, 2007)  run for 300 iterations within the OSTRICH optimization package  (Matott, 2017). A list of the calibrated nine parameters is given in Table  S3  in the Supplement.
All forcing data and geophysical basin characteristics used are provided in the common dataset (Sect. 2. 2). In total, 11 model parameters are calibrated against the KGE using the shuffled complex evolution (SCE) algorithm  (Duan et al. , 1992). The SCE was set up with two complexes and a maximum of 25 (internal) loops, and the algorithm converged on average after 1800 model evaluations. One calibration trial was performed for each basin. The list of the 11 calibrated parameters is provided in Table  S4  in the Supplement. The model variables for AET and SWE are vertically lumped and saved as model outputs to evaluate them against reference datasets. The proxy variable used for SSM is the storage content in the nonlinear storage tank (C) from which the runoff is derived (see the schematic diagram of HYMOD2 shown in Fig.
The model parameters calibrated can be found in Sect. S3. 7 and Table  S8  in the Supplement. -VIC-Raven. The Variable Infiltration Capacity (VIC) model  (Liang et al. , 1994; Hamman et al. , 2018 ) is a semi-distributed macroscale hydrologic model that has been widely applied in water and energy balance studies  (Newman et al. , 2017; Xia et al. , 2018; Yang et al. , 2019). In this study, VIC was forced with the gridded RDRS-v2 meteorological data, including precipitation, air temperature, atmospheric pressure, incoming shortwave radiation, incoming longwave radiation, vapor pressure, and wind speed at the 6 h time step and 10 km by 10 km spatial resolution. Other data, such as topography, soil, and land cover, are derived from those datasets, as described in Sect. 2. 2. VIC was used for runoff generation simulation at grid cell level, and the resultant daily gridded runoff and recharge fluxes were aggregated to subbasin level for the Raven routing, which then produced daily streamflow at the outlet of each catchment.
, 1998). For each of the six regions, the geophysical datasets (DEM, soil, and land use), as presented in Sect. 2. 2, were used to create subbasins which were further discretized into hydrologic response units. Subbasin-averaged daily forcing data (maximum and minimum temperature and precipitation), as presented in Sect. 2. 2, were then used simulate land-surface processes. The vertical flux (total runoff) produced at the subbasin spatial scale was then fed to a lake and river routing product  (Han et al. , 2020) , integrated into Raven modeling framework  (Craig et al. , 2020) , for subsequent routing. Note that the SWAT subbasins and routing product subbasins are not equivalent, and hence, aggregation to routing product subbasins was necessary, similar to the aggregation for gridded models.
A separate WATFLOOD model for each of the six regions (Fig. 1c ) was developed using the common geophysical datasets and hourly precipitation and temperature data at a regular rotated grid of about 0. 09 decimal degrees (≈ 10 km). Model grids were further discretized into grouped response units. Hourly gridded runoff and recharge to lower zone storage (LZS) were aggregated to the subbasin level for the Raven routing, which then produced hourly streamflow at the outlet of each catchment. Streamflow calibration and validation of the integrated WATFLOOD-Raven was carried out using the same methodology (different parameters) as described in the previous section for SWAT-Raven. We refer to Sect. S3. 10 and Table  S11  in the Supplement for further details of the 11 WATFLOOD and 6 Raven parameters used during model calibration.
g. ,  Haghnegahdar et al. , 2015). The MESH model setup here is using the Canadian Land Surface Scheme (CLASS) to produce grid-based model outputs, which are passed to the Raven routing to generate the streamflow at the locations used in this study. The model will hence be referred to as MESH-CLASS-Raven hereafter. The following meteorologic forcing variables from the RDRS-v2 dataset (Sect. 2. 2) are used to drive CLASS: incoming shortwave radiation, incoming longwave radiation, precipitation rate, air temperature, wind speed, barometric pressure, and specific humidity (all provided at or near surface level). The model simulates vertical energy and water fluxes for vegetation, soil, and snow. The MESH-CLASS-Raven grid is the same as the forcing data (≈ 10 km by 10 km) and produces outputs at an hourly time step.
MESH allows, for example, us to emulate the GEM-Hydro-Watroute model outside of ECCC informatics infrastructure. In this study, however, the hourly MESH-SVS gridded outputs of total runoff and drainage were provided to the Raven routing scheme, which is computationally less expensive than the Watroute routing scheme (more info in Sect. S3. 13 in the Supplement). The model will be referred to as MESH-SVS-Raven hereafter. The following meteorological forcing variables from the RDRS-v2 dataset (Sect. 2. 2) are used to drive the SVS LSS (same as that used for MESH-CLASS-Raven): incoming shortwave radiation, incoming longwave radiation, precipitation rate, air temperature, wind speed, barometric pressure, and specific humidity. SVS also needs geophysical fields such as soil texture, vegetation cover, slope, and drainage density.
All of these variables were derived from the common geophysical datasets used for this project (Sect. 2. 2), except for drainage density, which is derived from the National Hydro Network (NHN) dataset over Canada (Natural Resources Canada, 2020) and from the National Hydrography Dataset (NHD) dataset over the United States of America (USGS, 2021). In this study, the MESH-SVS-Raven grid is the same as the forcing data (around 10 km by 10 km), and SVS was run using a 10 min time step. The Raven routing was used with a 6 h time step. Routing results were aggregated to a daily time step for model performance evaluation. MESH-SVS-Raven was calibrated using a global calibration approach  (Gaborit et al. , 2015)  but applied for each of the six regions of the complete area of interest, therefore resulting in what is called in this study a regional calibration method.
The model will be referred to as GEM-Hydro-Watroute in this study. It was preferred to use its native gridded routing scheme of Watroute during this study, instead of the common Raven routing scheme used by other hydrologic models, in order to investigate the potential benefit of calibration for the experimental ECCC hydrologic forecasting system named the National Surface and River Prediction System (NSRPS;  Durnford et al. , 2021) , which relies on GEM-Hydro-Watroute. GEM-Surf forcings are the same as those used for MESH-SVS-Raven, except that GEM-Surf needs both wind speed and wind direction (available in RDRS-v2), and surface water temperature and ice cover fraction (taken from the analyses of the ECCC Regional Deterministic Prediction System). Regarding the geophysical fields needed by GEM-Surf, they consist of the same as those used for MESH-SVS-Raven, except that some additional fields are required, like surface roughness and elevation (taken or derived from the common geophysical datasets used in this project).
Moreover, the Watroute model used here mainly relies on Hy-droSHEDS 30 arcsec (≈ 1 km) flow direction data (Hy-droSHEDS, 2021), Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010; USGS, 2010), and on Climate Change Initiative Land Cover (CCI-LC) 2015 vegetation cover  (ESA, 2015). The surface component of GEM-Hydro-Watroute was employed here with the same resolution as the forcings (≈ 10 km) and with a 10 min time step, with outputs aggregated to hourly afterwards. The Watroute model used here has a 30 arcsec (≈ 1 km) spatial resolution and a dynamic time step comprised between 30 and 3600 s. Daily flow averages were used when computing model performances. GEM-Hydro-Watroute is computationally intensive; it was not optimized to perform long-term simulations but is rather oriented towards large-scale forecasting.
Model calibration and validation setup and datasets The models participating in this study were calibrated and validated based on streamflow observations made available by either Water Survey Canada (WSC) or the U. S. Geological Survey (USGS). The streamflow gauging stations selected for this project had to have at least a (reported) drainage area of 200 km 2 (to avoid too flashy watershed responses in the hydrographs) and less than 5 % of missing (daily) data between 2000 and 2017, which is the study period, i. e. , warmup period (January to December 2000), calibration period (January 2001 to December 2010), and validation period (January 2011 to December 2017). The low portion of missing data was to avoid not having enough data to properly estimate the calibration and validation performance of the models.
There are certainly more appropriate/advanced ways to do regionalization, but the main objective was to ensure that the models used the same approach. We wanted to avoid favoring a few models by testing more advanced approaches with them and then applying that method to other models which might have favored other strategies. The donor basin mapping is provided in Table S16 in the Supplement. The WSC and USGS streamflow data were downloaded by the project leads, the units were converted to be consistent in cubic meters per second (m 3 s -1 ), and then shared with the team in both CSV and NetCDF file formats to make sure all team members used the same version of the data for calibration. Note that at first only the data from the 141 calibration locations was shared. The validation locations and respective observations were only made available after the calibration of all models was finished to allow for a true blind (spatial and spatiotemporal) validation.
Model evaluation setup and datasets In order to assess the model performance beyond streamflow, a so-called model evaluation was performed. Model evaluation means that a model was not (explicitly) trained against any of the data used for evaluation. We distinguish model evaluation from the model validation of streamflow because, when building the models, although the modeling teams knew their model would be assessed based on streamflow validation, they did not know in advance that additional variables beyond streamflow would be assessed against other observations (this was only decided later in the study). It is not known if the models have been trained in previous studies against any of the observations we use for validation or model evaluation. Any model previously using such data to inform model structure or process formulations might have an advantage relative to another model whose structural development did not involve testing in this region.
Such kinds of model developments are, however, the nature of modeling in general and are assumed to determine the quality of a model that is to be assessed here exactly. The evaluation datasets were chosen based on the following criteria: (a) they had to be gridded to allow for a consistent evaluation approach across all additional variables, (b) the time step of the variables preferably had to be daily to match the time step of the primary streamflow variable, (c) they had to be available over the Great Lakes for a long time period between 2001 and 2017 (study period) to allow for a reliably long record in time and space, (d) they had to be open source to allow for the reproducibility of the analysis, (e) they had to be variables typically simulated in a wide range of hydrologic and land surface models, and (f) they had to have been at least partially derived on distributed observational datasets for the region (e.
g. , ground or satellite observations). There were two evaluation variables taken from the GLEAM (Global Land Evaporation Amsterdam Model) v3. 5b dataset  (Martens et al. , 2017). GLEAM is a set of algorithms using satellite observations of climatic and environmental variables in order to estimate the different components of land evaporation and surface and root zone soil moisture, potential evaporation, and evaporative stress conditions. The first variable used in this project is actual evapotranspiration (AET; variable E in the dataset; see Fig. 2b ) given in millimeters per day (mm d -1 ), and the second variable is surface soil moisture (SSM; variable SMsurf in the dataset; see Fig. 2c ) given in cubic meters per cubic meter (m 3 m -3 ) and valid for the first 10 cm of the soil profile.