At elevated temperatures, pMOS transistors show a considerable drift in fundamental device parameters such as the threshold voltage when a large negative bias is applied. This phenomenon, known as negative bias temperature instability, is regarded as one of the most important reliability concerns in highly scaled pMOS transistors. Modeling efforts date back to the reaction-diffusion (RD) model proposed by Jeppson and Svensson forty years ago which has been continuously refined since then. So far, the change in the interface state density predicted by the RD model is directly used to approximate the threshold voltage shift. Here we present a coupling of the RD model to the semiconductor equations which is required to go beyond that approximation and to study degradation during realistic device operating conditions. It is also shown that such a coupled treatment is required to accurately model the behavior during the measurement phase. In addition, the RD model is extended to improve the prediction both in the stress and the relaxation phase by accounting for trap-controlled transport of the released hydrogen species.