Gas-liquid liquid solid (GLLS) reaction systems are often encountered in manufacturing of fine and specialty chemicals. More often than not, such reaction systems involve multiple reactions, and selectivity toward the desired component always poses challenges. An adequate understanding of various parameters affecting GLLS reactor performance is essential to develop strategies for realizing desired selectivity. In this work, a comprehensive reaction engineering model for simulating four phase hydrogenation reactions has been developed. A generalized mixing cell based framework for a reaction system with four interacting phases (gas [G], aqueous [L], organic [L], and solid catalyst [S]) was developed. The model is written in a general way so as to specify one of the liquid phases as a continuous phase, and the other three phases are dispersed into it. In each cell, vapor space is included. The model includes the possibility of evaporation of solvent and internal condensation (in vapor space). The model can also be applied for a dead end (from a perspective of reacting gas) reactor. Model equations were solved using MATLAB. The equations and solution methodology were verified by comparing numerical solutions with available solutions of various limiting cases. A case of four phase hydrogenation of nitrobenzene to para amino phenol and aniline was considered to illustrate the application of the developed model. Key findings from the model were validated by comparing with laboratory scale experimental data. The model was then used to develop insights and guidelines for enhancing selectivity toward desired product. The developed model and presented results will be useful to develop general guidelines for design and optimization of GLLS reactors.