The sodium column abundance CNa in the mesosphere is one of the crucial parameters that determines the brightness of sodium laser guide stars (LGS). LGS are becoming standard elements of adaptive optics assisted instruments and more powerful lasers as well as systems to generate multiple LGS from one single laser are available now. An accurate knowledge of the possible LGS photon return flux is therefore important to properly dimension the laser power, for accurate performance estimation through adaptive optics simulations, and to design the LGS wavefront sensors. We have been using the telemetry data of the 4-LGS Adaptive Optics Facility at Paranal available since mid-2017 to derive CNa. The automatic logging provides us with an almost unique data set spanning already more than five years and counting. We derive the statistics and variations of CNa on different time scales: hours, days, months, seasons, years, or even decades. With these data, we developed a semi-empirical model of CNa at Cerro Paranal. Moreover, we propose here a global model of CNa, for every location and any time. To do so, we have gathered additional data from various sources and locations, coming from other observatories, Lidars, Earth observation satellites, and Earth atmosphere models. We will present the status of this modelling attempt, showing the challenges and the careful verifications we applied to our data analysis to avoid biases, and give a glimpse at a public domain numerical prediction tool.