The passive scalar spreading of fluids with laminar Prandtl or Schmidt number, , Pr Sc , equal to 1 in turbulent rectangular submerged free jets is analyzed by means of numerical simulation and theoretical analysis in the Reynolds number range 5000-40,000. The numerical investigation is carried out by means of a three-dimensional (3D) Large Eddy Simulation (LES) approach with the dynamic Smagorinsky model. A new mathematical model allows to obtain a simplified description of the passive scalar spreading in the largest area of the flow field, the Fully Developed Region (FDR). The present three-dimensional (3D) investigation shows that the passive scalar spreading follows a self-similarity law in the Fully Developed Region (FDR), as well as in the mean Undisturbed Region of Flow (URF) and in the Potential Core Region (PCR), similarly to what found in the Near Field Region (NFR) of rectangular submerged free jets, investigated with a two-dimensional (2D) approach. The turbulent Prandtl or Schmidt number is evaluated numerically and is found to be inversely proportional to the mean velocity gradient in the PCR. The present 3D numerical results show that the turbulent Prandtl or Schmidt number is zero in most part of the mean URF, and PCR, while it assumes different values outside. In the FDR the turbulent Prandtl or Schmidt number is constant and approximately equal to 0.7, in agreement with the literature, showing that turbulence affects momentum and passive scalar in a different way.