Advancements in computing resources continually motivate increases in the complexity and spatial resolution of global atmospheric simulations. However, the impact of the numerical methods used for coupling different processes on the fidelity of the simulations is often neglected. While non-physical model behavior and undesirable sensitivities can sometimes be alleviated by relatively simple code changes based on atmosphere modelers’ intuition, a more systematic approach to the selection of numerical methods is needed to faithfully represent the complex interactions among atmospheric processes in the real world. Our team, funded by DOE’s SciDAC program, has been working on this challenge through collaborations between Earth scientists and applied mathematicians. This presentation will demonstrate that, by combining physical understanding with mathematical error analysis, one can obtain insights into the pros and cons of numerical coupling schemes used in current models, as well as design alternative schemes aiming at good trade-offs between numerical accuracy, computational efficiency, and software simplicity. We will report on our progress using examples related to clouds and aerosols in the E3SM atmosphere model.