This work develops an integrated method for characterizing the spatiotemporal variability and distribution function of extreme total water level events in a transitional coastal environment forced by multiple simultaneous climatic agents. The approach incorporates the reconstruction of the total and isolated components of the water level, as well as the nonlinear interaction term, via the hybrid downscaling technique. It also analyzes the contribution of both non-extreme and extreme values of each component to the magnitude and variability of the extreme total water level events and their relations of concomitance to characterize the combinations that cause the extreme events. In addition, the method assesses the probability of flooding events along the transitional coastal environment. The outcomes of the method constitute a useful tool for dividing the transitional environment into homogeneous zones based on the probability of flooding, the dominance of the agents and the relationships between water level components. The application of the proposed method to the Guadalete estuary (South West, Spain) enables the identification of three homogeneous sections: the outer section, where there is no significant variability of the extreme total water level events, and they are explained only by the extremes of the astronomical tide; the inner section, with the highest variability and the domination of the extremes of the river discharge over the rest of the agents, and the middle estuary, where a similar contribution between the extreme values of the river discharge and the extreme and non-extreme values of the astronomical tide is found. These zonation maps are key elements to support decision-making for ecosystem management and risk analysis and facilitate the identification of vulnerable areas, quantification of the flooding frequency and identification of the agents responsible for flooding.