The co-occurrence of (not necessarily extreme) precipitation and surge can lead to extreme inland water levels in coastal areas. In a previous work the positive dependence between the two meteorological drivers was demonstrated in a managed water system in the Netherlands by empirically investigating an 800-year time series of water levels, which were simulated via a physical-based hydrological model driven by a regional climate model large ensemble. In this study, we present an impact-focused multivariate statistical framework to model the dependence between these flooding drivers and the resulting return periods of inland water levels. This framework is applied to the same managed water system using the aforementioned large ensemble. Composite analysis is used to guide the selection of suitable predictors and to obtain an impact function that optimally describes the relationship between high inland water levels (the impact) and the explanatory predictors. This is complex due to the high degree of human management affecting the dynamics of the water level. Training the impact function with subsets of data uniformly distributed along the range of water levels plays a major role in obtaining an unbiased performance. The dependence structure between the defined predictors is modelled using two- and three-dimensional copulas. These are used to generate paired synthetic precipitation and surge events, transformed into inland water levels via the impact function. The compounding effects of surge and precipitation and the return water level estimates fairly well reproduce the earlier results from the empirical analysis of the same regional climate model ensemble. Regarding the return levels, this is quantified by a root-mean-square deviation of 0.02gm. The proposed framework is able to produce robust estimates of compound extreme water levels for a highly managed hydrological system. Even though the framework has only been applied and validated in one study area, it shows great potential to be transferred to other areas. In addition, we present a unique assessment of the uncertainty when using only 50 years of data (what is typically available from observations). Training the impact function with short records leads to a general underestimation of the return levels as water level extremes are not well sampled. Also, the marginal distributions of the 50-year time series of the surge show high variability. Moreover, compounding effects tend to be underestimated when using 50-year slices to estimate the dependence pattern between predictors. Overall, the internal variability of the climate system is identified as a major source of uncertainty in the multivariate statistical model.