Multilevel models with nested errors are widely used in poverty estimation. An important application in this context is estimating the distribution of poverty as defined by the distribution of income within a set of domains that cover the population of interest. Since unit-level values of income are usually heteroskedastic, the standard homoskedasticity assumptions implicit in popular multilevel models may not be appropriate and can lead to bias, particularly when used to estimate domain-specific income distributions. This article addresses this problem when the income values in the population of interest can be characterized by a two-level mixed linear model with independent and identically distributed domain effects and with independent but not identically distributed individual effects. Estimation of poverty indicators that are functionals of domain-level income distributions is also addressed, and a nonparametric bootstrap procedure is used to estimate mean squared errors and confidence intervals. The proposed methodology is compared with the well-known World Bank poverty mapping methodology for this situation, using model-based simulation experiments as well as an empirical study based on Bangladesh poverty data.