The Groningen gas field is one of the largest gas fields in Europe. The continuous gas extraction led to an induced seismic activity in the area. In order to monitor the seismic activity and study the gas field many permanent and temporary seismic arrays were deployed. In particular, the extraction of the shear wave velocity model is crucial in seismic hazard assessment. Local S-wave velocity-depth profiles allow the estimation of a potential amplification due to soft sediments. Ambient seismic noise tomography is an interesting alternative to traditional methods that were used in modelling the S-wave velocity. The ambient noise-field consists mostly of surface waves, which are sensitive to the S-wave and if inverted, they reveal the corresponding S-wave structures. In this study, we present results of a depth inversion of surface waves obtained from cross-correlation of 1 month of ambient noise data from four flexible networks located in Groningen area. Each block consisted of 400 3C-stations. We compute group velocity maps of Rayleigh and Love waves using a straight-ray surface wave tomography. We also extract clear higher modes of Love and Rayleigh waves. The S-wave velocity model is obtained with a joint inversion of Love and Rayleigh waves using the Neighbourhood Algorithm. In order to improve the depth inversion, we use the mean phase velocity curves and the higher modes of Rayleigh and Love waves. Moreover, we use the depth of the base of the North Sea formation as a hard constraint. This information provides an additional constraint for depth inversion, which reduces the S-wave velocity uncertainties. The final S-wave velocity models reflect the geological structures up to 1 km depth and in perspective can be used in seismic risk modelling.