Accurate characterization and modeling of subsurface flow in multi-fluid soil systems require development of a rapid, reproducible experimental method that yields the information necessary to determine the parameters of the capillary pressure-saturation and permeability functions. Previous work has demonstrated that parameter optimization using inverse modeling is a powerful approach to indirectly determine these constitutive relationships for air-water systems. We consider extension of the inverse parameter estimation method to the modified multi-step outflow method for two-fluid (air-water, air-oil and oil-water) flow systems. The wellposedness of the proposed parameter estimation problem is evaluated by its accuracy, uniqueness and parameter uncertainty. Seven different parametric models for the capillary pressuresaturation and permeability functions were tested in their ability to fit the multi-step outflow experimental data; these included the van Genuchten-Mualem (VGM) model, van Genuchten-Burdine (VGB), Brooks-Corey-Mualem (BCM), Brooks-Corey-Burdine (BCB), Brutsaert-Burdine (BRB), Gardner-Mualem (GDM) and Lognormal Distribution-Mualem (LNM) models. The VGB, BCM and BCB models fitted the multi-step outflow data poorly, and resulted in relatively large values for the rootmean-squared residuals. It was concluded that the remaining VGM, LNM, BRB and GDM models successfully characterized the multi-step experimental data for two-fluid flow systems. Although having one parameter less than the other models, the GDM model's performance was excellent. Finally, we conclude that optimized capillary pressure-saturation functions for the oil-water and air-oil could be predicted from the respective air-water function using interfacial tension ratios.