Summary Poverty maps at local level might be misleading when based on direct (or area-specific) estimators obtained from a survey that does not cover adequately all the local areas of interest. In this case, small area estimation procedures based on assuming common models for all the areas typically provide much more reliable poverty estimates. These models include area effects to account for the unexplained between-area heterogeneity. When poverty figures are sought at two different aggregation levels, domains and subdomains, it is reasonable to assume a twofold nested error model including random effects explaining the heterogeneity at the two levels of aggregation. The paper introduces the empirical best (EB) method for poverty mapping or, more generally, for estimation of additive parameters in small areas, under a twofold model. Under this model, analytical expressions for the EB estimators of poverty incidences and gaps in domains or subdomains are given. For more complex additive parameters, a Monte Carlo algorithm is used to approximate the EB estimators. The EB estimates obtained of the totals for all the subdomains in a given domain add up to the EB estimate of the domain total. We develop a bootstrap estimator of the mean-squared error of EB estimators and study the effect on the mean-squared error of a misspecification of the area effects. In simulations, we compare the estimators obtained under the twofold model with those obtained under models with only domain effects or only subdomain effects, when all subdomains are sampled or when there are unsampled subdomains. The methodology is applied to poverty mapping in counties of the Spanish region of Valencia by gender. Results show great variation in the poverty incidence and gap across the counties from this region, with more counties affected by extreme poverty when restricting ourselves to women.