This study aims to investigate how uncertainty from model parameters and physics parameterizations propagates through a regional climate model of the Great Lakes Region (GLR) that couples the atmosphere with lake hydrodynamics. We introduce a generalized framework that analyzes the full spatiotemporal variation of the quantities of interest (QoI) using a relatively small ensemble of the high-dimensional mixed parameter/parameterization (continuous/discrete) space. In the current application, there are 9 such dimensions: 4 parameterizations of the atmosphere, 2 parameters and 2 parameterizations of the lake, and 1 parameterization of the land surface, and 18 model ensembles using quasi-random sequence sampling. Principal component analysis is performed on the full spatiotemporal distribution of QoIs: surface air and lake temperatures and wind/current speeds, and surrogate representations of the mode of variations are generated using Polynomial Chaos Expansions (PCE) and Neural Network (NN) models. The study highlights how the parsimonious and analytical representation by PCEs can be beneficial over NN for analyzing large spatiotemporal scales. We show how parameterizations of the atmosphere affect lake temperatures (see included figure) and current speeds and vice-versa, and interactions between the parameters/parameterizations. We see different behavior between the shallower and deeper lakes and remaining structural error of the temperature distribution in the deep lakes during the early summer warming phase.