Abstract: We study the production of the four-lepton final state l+l−l+l−, predominantly produced by a pair of electroweak Z bosons, ZZ. Using the LoopSim method, we merge NLO QCD results for ZZ and ZZ+jet and obtain approximate NNLO predictions for ZZ production. The exact gluon-fusion loop-squared contribution to the ZZ process is also included. On top of that, we add to our merged sample the gluon-fusion ZZ+jet contributions from the gluon-gluon channel, which is formally of N^3LO and provides approximate results at NLO for the gluon-fusion mechanism. The predictions are obtained with the VBFNLO package and include the leptonic decays of the Z bosons with all off-shell and spin-correlation effects, as well as virtual photon contributions. We compare our predictions with existing results for the total inclusive cross section at NNLO and find a very good agreement. Then, we present results for differential distributions for two experimental setups, one used in searches for anomalous triple gauge boson couplings, the other in Higgs analyses in the four charged-lepton final state channel. We find that the approximate NNLO corrections are large, reaching up to 20% at high transverse momentum of the Z boson or the leading lepton, and are not covered by the NLO scale uncertainties. Distributions of the four-lepton invariant mass are, however, stable with respect to QCD corrections at this order.