Introduction: In this study, the two-dimensional equations of energy, mass, and momentum were solved using Comsol Multiphysics 5.2 software for pasteurizing date syrup in polyethylene terephthalate (PET) packaging. The importance of date syrup has been quite studied before in terms of its valuable nutritional values (especially heat-sensible vitamins). Research has shown that the appropriate heating process is required to maintain these values. Modeling (especially numerical modeling) is a robust method to predict the product temperature profile at every sample point during the process. It accurately predicts the slowest heating zone (SHZ), which should be reached at the PASTEURIZATION temperature and remain quiet enough at that temperature. In this way, not only are resistant microorganisms deactivated and killed but the product's nutritional values are also maintained. Comsol Multiphysics was applied in this research to predict the accurate position of the slowest heating zone, the dominant heat transfer method, and fluid velocity.Materials and methods: Thermophysical properties are necessary for heat transfer models. We estimated thermal conductivity (k) using the Krischer model. This model estimates a more logical thermal conductivity because both series and parallel models are incorporated into this model. The ratio of parallel and series models was assumed to be the same (f = 0.5) in our research. Specific heat was modeled using a suit model based on the mass fraction of components and their specific heat by temperature through a parallel model. Density is measured using a defined volume container. For this purpose, a pycnometer was used, and its volume was measured using distilled water filled in. The weight of syrup divided into its volume gave the density. The surface heat transfer coefficient (h) was determined by unsteady temperature measurements. This method measured the temperature of an aluminum container exactly at the same size as the PET, and h was obtained using the slope of ln(T ±Tal) diagram. The container's geometry should be drawn at the first modeling stage. In our case, the container was a hard part as the container did not have a specific geometric shape and had a thickness with thermal resistance. Second, the physics based on appropriate equations were selected, followed by defining initial and boundary conditions. The density, thermal conductivity, and surface heat transfer values were 1376 kg/m3, 0.4 W/m℃, and 43 W/m2℃, respectively, and specific heat was temperature-dependent. Next, meshes were defined. In the present study, the number of meshes used in the model was 6532 triangular elements, 405 of which were in boundaries. There should be logical assumptions to be able to model a process. In our study, the assumptions of Ghani were used. Finally, the model was run. In the experimental part of this research, heating was performed at 70℃ for 3000 s. Then, the temperatures of different container parts and syrup were collected using a K-type thermocouple and a data logger. Finally, the collected data at different parts of the package were used to verify the model.Results and discussion: The comparison criteria between the predicted and experimental figures used to evaluate the goodness of fit (GoF), namely the correlation coefficient (r =0.999), indicated that the model was valid and we could benefit from the model results. The results showed that the cold point migrated towards the top of the container because of the high product viscosity and the big can geometry. In fact, the dominant way of heat transfer was conduction. The explanation is that when heat transfer is molecule-by-molecule, the position of the cold point migrates toward the top of the container as at the top boundaries. Therefore, vacuum condition at headspace is considered. The accurate position of the slowest heating zone was at r = 9.5 cm and z = 0. The fluid velocity was maximum near the can wall and at a lower height and minimum at the beginning and the end of the process due to less temperature gradient.The mean velocity value in this situation was 1.562 ×10 -7 m/s. The plots of fluid velocity versus container radius and height showed that by increasing the height, fluid velocity rose because the fluid warmed up and flew toward the top of the can. The velocity in the interior radius after about 50 min of the heating part was more than the wall vicinity because the warmed fluid migrated toward the top of the can would be immobilized, and heat exchange with the cold parts would occur there. As a result, the fluid with a lower temperature would return to the bottom. In this process, after about 50 min when the fluid around the wall reached the environment temperature, there would be a temperature gradient near the core with colder fluid and the wall, which resulted in higher fluid velocity. The results of fluid velocity during the heating time show that this factor would decrease by increasing heating time due to the reduction in a temperature gradient. Obviously, the minimum velocity would be at the beginning and the end of the heating process when the temperature gradient was minimum. Different velocity plots demonstrated that at the beginning times of the heating process, the maximum velocity was near the wall (r = 4.5 cm). In addition, at the end of the process, it was at the interior parts near the center (r = 0.5 cm) regarding the difference in the magnitude of the temperature gradient at different parts of the fluid. We can also compare heat transfer by conducting or combining conduction-convection in our model. The model showed that the maximum heat flux is 320 W/m2 for conduction and 30 W/m2 for the convection part. These figures indicate the importance of conduction in the date syrup sample is more than ten times convection. The magnitude of heat flux in both conditions after around 35 min and by temperature stabilization at 70℃ reduced dramatically. Conclusion: The simulation showed that the required time for PASTEURIZATION was 35 min, and the cold point reached the autoclave temperature.