The reliability of fine pitch solder joints in advanced microelectronic packaging is often compromised by the growth of intermetallic compounds (IMC) along Cu pillar sidewalls, accompanied by solder voiding. Modeling the void growth requires accounting for surface diffusion that is the primary contributor to the observed phenomenon. Relatively few simulations in the literature account for surface diffusion. In this work, a reaction-diffusion framework incorporating both surface and bulk diffusion mechanisms is used to model IMC growth and void evolution in solder microbumps. The governing equations are derived using rigorous first-principles continuum mechanics that describe the balance of mass, linear momentum, and thermodynamics on an interface. A governing condition for interface velocity is then derived considering gradients in concentration, temperature, and stress. The void growth is simulated using the previously developed Enriched Isogeometric Analysis (EIGA) method, which allows explicit interface tracking through parametric splines, facilitating the direct implementation of the governing equations. There is no need for additional numerical stabilization or complex modified governing equations which are often necessary in phase field and level set methods. The boundary conditions at the interface are applied exactly in EIGA. The simulations provide insights into the thermal aging on void evolution rates and assess the reliability of solder joints across different bump sizes and aging conditions.