The fitting algorithm in gstat is only one of several possibilities, in particular if using weighted least squares then there are several different choices for the weightings. It is not really a question of how well the variogram model fits the variogram cloud but rather how well the variogram model fits the data (the data and data locations determine the variogram cloud but not conversely). Cross validation is one way to judge how well the variogram model fits the data and there are multiple statistics for that; mean error, normalized mean square error, fraction of normalized errors exceeding 2.5 (in absolute value). Since kriging is an unbiased estimator the expected values of the mean error is zero, likewise the expected value of the normalized mean square error is one, the probability of the the fraction of normalized errors exceeding 2.5 is less than .05 (using Chebyschev's inequality). It is also useful to examine the histogram of the errors and of the normalized errors as well as a coded plot of the data locations (coded by the normalized errors). No one of these is an absolute determinant. AIC and BIC are really based on different assumptions and not really applicable. Fitting a variogram is not like fitting a linear regression model