Quenching of heated surfaces through impinging liquid jets is of great im-portance for numerous applications like steel processing, nuclear power plants, automobile industries, etc. Therefore, computational modelling of the surface quenching through circular water jets impinging normally onto a heated flat surface has vital importance in order to reveal the physics of the quenching process. At first, a numerical model was developed for single jet impingement process. A conjugate heat transfer problem was solved implying consideration of both regions, one occupied by fluid (multi-phase flow consisting of water, vapor and ambient air) and one accommodating the solid surface within the same solution domain. Numerical simulations were performed in a range of relevant operating param-eters: jet velocities (2.5, 5, 7.5 and 10 m/s), water sub-cooling (75 K) and wall-superheat (650 K - 800 K) corresponding closely to those encountered in the industrial water jet cooling banks. Due to the high initial temperature of the surface, the boiling process exhibits strong spatial and temporal fluctuations. Its effect on the boundary layer profiles at the stagnation region at different time intervals are analyzed. The analysis reveals a highly distorted field of both mean flow and turbulence quantities. It represents an im-portant outcome, also with respect to appropriate model improvements. The different boiling characteristics are envisaged in detail to increase the level of understanding of the phenomena. The influence of the turbulent kinetic energy investigated at the boiling front as well as the jet-acceleration region has been studied. The physically relevant results are obtained and analyzed along with reference database provided experimentally by Karwa (2012, ‘Experimental Study of Water Jet Impingement Cooling of Hot Steel Plates,' Dissertation, FG TTD, TU Darmstadt). The intensive quenching process is consistent with the high rate of sub-cooling and high jet velocities. The surface temperature predicted by quenching model within the impingement region and the subsequent wall-jet region agrees reasonably well with the measurements, the outcome being particularly valid at higher jet velocities. However, a steep temperature gradient at the position corresponding to boiling threshold has not been captured under the condition of the numerical grid adopted. On the other hand, a reasonably good prediction of the wetting front propagation phenomena advocates the future development of the model. The high-intensity back motion of the vapor phase in the stagnation region at the earlier times of the water jet impingement can induce an appropriately high turbulence level, which could be accounted well by the turbulence model applied. The second part of the present work deals with multiple liquid jet impingement. When the multiple jets impact onto the heated surface, the heat flux is extracted from the surface by the mass flow rate. The heat flux is dependent on the several flow conditions and configurations of the nozzle array system. Therefore, one needs to study the nozzle array configuration along with several flow parameters for the better design of the cooling header system. Accordingly, the hydrodynamics of the multiple jets has been studied computationally realizing the need for optimum configuration of the nozzle array. The effects of the mass flow rate, target plate width and the turbulence produced due to the impingement were studied. Afterward, an analytical model is proposed for the quenching of the multiple jets system. It has been realized that, when jet impinges onto the surface at very high initial temperature, the film boiling may play a role in the heat transfer mechanism. Therefore, development has been made in the film-boiling model considering the effect of turbulence at the liquid jet stagnation region at the Leidenfrost point. The Leidenfrost point is the minimum temperature at which the film boiling can sustain. However, the vapor-liquid interface has the dynamic character; it oscillates with high frequency and causes the additional momentum diffusivity. Therefore, the need for introducing the effect of associated turbulence has been felt. The length and velocity scale of the turbulent structure has been approximated by assuming homogeneous turbulence. The new model for the heat flux and wall superheat yielded results agreeing well with published experimental results.