Observations from the Interface Region Imaging Spectrograph (IRIS) often reveal significantly broadened and non-reversed profiles of the Mg II h, k and triplet lines at flare ribbons. To understand the formation of these optically thick Mg II lines, we perform plane parallel radiative hydrodynamics modeling with the RADYN code, and then recalculate the Mg II line profiles from RADYN atmosphere snapshots using the radiative transfer code RH. We find that the current RH code significantly underestimates the Mg II h & k Stark widths. By implementing semi-classical perturbation approximation results of quadratic Stark broadening from the STARK-B database in the RH code, the Stark broadenings are found to be one order of magnitude larger than those calculated from the current RH code. However, the improved Stark widths are still too small, and another factor of 30 has to be multiplied to reproduce the significantly broadened lines and adjacent continuum seen in observations. Non-thermal electrons, magnetic fields, three-dimensional effects or electron density effect may account for this factor. Without modifying the RADYN atmosphere, we have also reproduced non-reversed Mg II h & k profiles, which appear when the electron beam energy flux is decreasing. These profiles are formed at an electron density of ∼ 8 × 10 14 cm −3 and a temperature of ∼ 1.4 × 10 4 K, where the source function slightly deviates from the Planck function. Our investigation also demonstrates that at flare ribbons the triplet lines are formed in the upper chromosphere, close to the formation heights of the h & k lines.