Thanks Andreas for providing us a beautiful and clean solution. I am providing to the community the final version of the code:
function θ_CALIBRATION(∑T, N_iT, N_iZ, θ, θcalibr, discret)
# For every depth we have measurements
X = θcalibr.∑T[1:θcalibr.Nt] / (60.0 * 60.0)
# X = θcalibr.Date[1:θcalibr.Nt]
X2 = ∑T[1:N_iT] / (60.0 * 60.0)
Style = ["smooth, red, very thick", "smooth, orange, very thick", "smooth, brown, very thick", "smooth, gray, very thick", "smooth, black, very thick"]
Plot_1 = map(axes(θcalibr.Data,2)) do iDepth
Label = "Neutron=" * string(discret.Znode[iDepth]) * "mm"
PGFPlots.Plots.Linear(X, θcalibr.Data[1:θcalibr.Nt, iDepth], legendentry=Label, style=Style[iDepth], mark="none")
end # loop
Plot_θcalibr = GroupPlot(1, 1, groupStyle = "horizontal sep = 2.5cm, vertical sep = 1.5cm")
push!(Plot_θcalibr, Axis(Plot_1, style="width=20cm, height=10cm", xlabel=L"$\ Time \ [Hour]$", ylabel=L"$\theta \ [cm^3 \ cm^{-3}]$", xmin=0.0, title=path.SiteName_Hypix, legendStyle ="{at={(0.2,-0.3)}, anchor=south west, legend columns=3}"))
PGFPlots.save(path.Hypix_θcalibr, Plot_θcalibr)
return
end # function: θ_CALIBRATION
