It is becoming widely accepted that catalytic promiscuity, i.e., the ability of a single enzyme to catalyze the turnover of multiple, chemically distinct substrates, plays a key role in the evolution of new enzyme functions. In this context, the members of the alkaline phosphatase superfamily have been extensively studied as model systems in order to understand the phenomenon of enzyme multifunctionality. In the present work, we model the selectivity of two multiply promiscuous members of this superfamily, namely the phosphonate monoester hydrolases from Burkholderia caryophylli and Rhizobium leguminosarum. We have performed extensive simulations of the enzymatic reaction of both wild-type enzymes and several experimentally characterized mutants. Our computational models are in agreement with key experimental observables, such as the observed activities of the wild-type enzymes, qualitative interpretations of experimental pH-rate profiles, and activity trends among several active site mutants. In all cases the substrates of interest bind to the enzyme in similar conformations, with largely unperturbed transition states from their corresponding analogues in aqueous solution. Examination of transition-state geometries and the contribution of individual residues to the calculated activation barriers suggest that the broad promiscuity of these enzymes arises from cooperative electrostatic interactions in the active site, allowing each enzyme to adapt to the electrostatic needs of different substrates. By comparing the structural and electrostatic features of several alkaline phosphatases, we suggest that this phenomenon is a generalized feature driving selectivity and promiscuity within this superfamily and can be in turn used for artificial enzyme design.