Tephra fallout hazard assessment is commonly undertaken with the development of probabilistic maps that rely on numerical models. Among the steps for map production, the definition of input parameters of the model (including atmospheric conditions), the physical approximations of the numerical simulations, and the probabilities of occurrence of different eruption types in specific time frames are among the most critical sources of uncertainty. In this paper, we present a tephra fallout hazard assessment study for two volcanoes (Cotopaxi and Guagua Pichincha) in Ecuador. We utilize the coupled PLUME-MoM/HYSPLIT models, and we develop a procedure for uncertainty quantification where: i) we quantify the uncertainty on eruptive source parameters and eruption type occurrence through expert elicitation; ii) we implement a new procedure for correlations between the different parameters, and iii) we quantify the uncertainty of the numerical model by testing it with past eruptions and by deriving coefficients of mean model overestimation/underestimation. Probability maps of exceedance, given a deposit thickness threshold, and thickness maps, given a probability of exceedence, are produced for eruption of sub-Plinian and Plinian types, which are then merged into single maps concerning the next eruption. These are described according to the uncertainty distribution of eruption type occurrence probabilities, in terms of their 5^th percentile, mean and 95^th percentile values. We finally present hazard curves describing exceeding probabilities in 10 sensitive sites within the city of Quito. Additional information includes the areal extent and the people potentially affected by different isolines of tephra accumulation.