In the present research, a numerical modeling approach of the initial stage of consolidation during spark plasma sintering on the microscopic scale is presented. The solution of a fully coupled thermo-electro-mechanical problem also accounting for grain boundary and surface diffusion is found by using a staggered way. The finite-element method is applied for solving the thermo-electro-mechanical problem while the finite-difference method is applied for the diffusion problem. A Lagrange-based non-linear formulation is used to deal with the detailed description of plastic and creep strain accumulation. The numerical model is developed for simulating the structural evolution of the involved particles during sintering of powder compacts taking into account both the free surface diffusion of the particles and the grain boundary diffusion at interparticle contact areas. The numerical results obtained by using the two-particle model—as a representative volume element of the powder—are compared with experimental results for the densification of a copper powder compact. The numerical and experimental results are in excellent agreement.