Weak vector boson scattering (VBS) at the LHC provides an excellent source of information on the structure of quartic gauge couplings and possible effects of physics beyond the SM in electroweak symmetry breaking. Parameterizing deviations from the SM within an effective field theory at tree level, the dimension-8 operators, which are needed for sufficiently general modeling, lead to unphysical enhancements of cross sections within the accessible energy range of the LHC. Preservation of unitarity limits is needed for phenomenological studies of the VVjj events which signify VBS. Here we develop a numerical unitarization scheme for the full off-shell VBS processes and apply it to same-sign W scattering, i.e. processes like qq→qqW+W+. The scheme is implemented within the Monte Carlo program VBFNLO, including leptonic decay of the weak bosons and NLO QCD corrections. Distributions differentiating between higher dimensional operators are discussed.