We evaluate the phase-space integrals that arise in double real emission diagrams for semi-inclusive deep-inelastic scattering at next-to-next-to-leading order in QCD. Utilizing the reverse unitarity technique, we convert these integrals into loop integrals, allowing us to employ integration-by-parts identities and reduce them to a set of master integrals. The master integrals are then solved using the method of differential equations and expressed in terms of Goncharov polylogarithms. By examining the series expansion in the dimensional regulator, we discover additional relations among some of the master integrals. As an alternative approach, we solve the master integrals by decomposing them into angular and radial components. The angular parts are evaluated using Mellin-Barnes representation, while special attention is given to the singular structures of the radial integrals to handle them accurately. Here the results are provided in terms of one-fold integrals over classical polylogarithms. This approach provides a clearer understanding of the origin of soft and collinear singularities.