Modeling and dynamic analysis of liquid extraction columns are essential for the design, control strategies and understanding of column behavior during start up and shutdown. Because of the discrete character of the dispersed phase, the population balance modeling framework is needed. Due to the mathematical complexity of the full population balance model, it is still not feasible for dynamic and online control purposes. In this work, a reduced mathematical model is developed by applying the concept of the primary and secondary particle method (Attarakih et al., 2009b, Solution of the population balance equation using the one primary and one secondary particle method (OPOSPM), Computer Aided Chemical Engineering, vol. 26, pp. 1333–1338). The method is extended to solve the nonhomogenous bivariate population balance equation, which describes the coupled hydrodynamics and mass transfer in an RDC extraction column. The model uses only one primary and one secondary particles, which can be considered as Lagrangian fluid particles carrying information about the distribution as it evolves in space and time. This information includes averaged quantities such as total number, volume and solute concentrations, which are tracked directly through a system of coupled hyperbolic conservation laws with nonlinear source terms. The model describes droplet breakage, coalescence and interphase solute transfer. Rigorous hyperbolic analysis of OPOSPM uncovered the existence of four waves traveling along the column height. Three of these are contact waves, which carry volume and solute concentration information. The dynamic analysis in this paper reveals that the dominant time constant is due to solute concentration in the continuous phase. On the other hand, the response of the dispersed phase mean properties is relatively faster than the solute concentration in the continuous phase. Special shock capturing method based on the upwind scheme with flux vector splitting is used, with explicit wave speeds, as a time–space solver. The model shows a good match between analytical and numerical results for special steady state and dynamic cases as well as the published steady state experimental data.