Results are presented for the nucleon isovector electromagnetic form factors using 11 ensembles generated by the MILC collaboration using the 2+1+1-flavors HISQ action. They span 4 lattice spacings $a \sim$ 0.06, 0.09, 0.12 and 0.15~fm and 3 values of $M_\pi \sim 135, 225$ and 315 MeV. High-statistics estimates are used to perform a simultaneous extrapolation in the lattice spacing, lattice volume and light-quark masses. The $Q^2$ dependence over the range 0.05-1.4 ${\rm GeV}^2$ is investigated using both the $z$-expansion and the dipole form. Final $z$-expansion estimates for the isovector r.m.s. radius are $r_E = 0.769(27)(30)$ fm $r_M = 0.671(48)(76)$ fm and $\mu^{p-n} = 3.939(86)(138)$ Bohr magneton. The first error is the combined uncertainty from the leading-order analysis, and the second is an estimate of the additional uncertainty due to using the leading order chiral-continuum-finite-volume fits. The dipole estimates, $r_E = 0.765(11)(8)$ fm, $r_M = 0.704(21)(29)$ fm and $\mu^{p-n} = 3.975(84)(125)$, are consistent with those from the $z$-expansion but with smaller errors. Our analysis highlights three points. First, all data from the eleven ensembles and existing lattice data on, or close to, physical mass ensembles from other collaborations collapses more clearly onto a single curve when plotted versus $Q^2/M_N^2$ as compared to $Q^2$ with the scale set by quantities other than $M_N$. The difference between these two analyses is indicative of discretization errors, some of which presumably cancel when the data are plotted versus $Q^2/M_N^2$. Second, the size of the remaining deviation of this common curve from the Kelly curve is small and can be accounted for by statistical and possible systematic uncertainties. Third, to improve lattice estimates, high statistics data for $Q^2 < 0.1$ ${\rm GeV}^2$ are needed.