The joint probability analysis of river water temperature (RWT) and low flow (LF) characteristics is essential as their combined effect can negatively affect aquatic species, e.g., ectotherm fish. Traditional multivariate models may not be as effective as copula-based methodologies. This study introduces a new multivariate approach, the nonparametric copula density framework, free from any distribution assumption in their univariate margins and copula joint density. The proposed framework utilized RWT and LF datasets collected at five different river stations in Switzerland. The study evaluates a nonparametric Gaussian kernel with six bandwidth selectors to model marginal densities. It employs nonparametric-based Beta kernel density, Bernstein estimator, and Transformation kernel estimator to approximate copula density with nonparametric and parametric margins. The performance of some parametric copula densities was also compared. The most justifiable models were employed to estimate bivariate joint exceedance probabilities and return periods (RPs). The Beta kernel copula with Gaussian kernel margins outperformed other models for most stations; Bernstein and Transformation copula with Gaussian kernel margins were better for only one station each. The univariate RPs (RWT or LF) are lower than the AND-joint but higher than OR joint case. As the percentile value of LF events (serve as a conditioning variable) increases, the bivariate joint RPs of RWT also increase. Higher values in RWT events result in higher RPs than lower values at the fixed percentile value of LF. All such estimated risk statistics are beneficial to analyze their mutual risk in aquatic habitats and freshwater ecosystems.