I can see calculating both as a quality assurance check. But if they differ by any amount (even a small amount), then why would you choose one over the other? My point is that your objective (prior to looking at the fit) should be used to decide which number to use. Otherwise you just have two different estimates.

Using the 202x202 dataset I use the following:

data = Import["Z2c.dat"];
z2 = Flatten[Table[{x, y, data[[x, y]]}, {x, 202}, {y, 202}], 1];
ListPlot3D[z2, PlotRange -> {All, All, {0, 100}}, BoxRatios -> {1, 1, 1}, AxesLabel -> Map[Style[#, Bold, Large] &, {"x", "y", "z"}]]

to obtain

I've truncated the z-axes to 100 to show the "noise" surrounding the "signal" in the middle. As such this is not about estimating a bivariate probability function but rather a function that has a surface similar to a bivariate normal. Because of the noise surrounding the signal I've added a constant to the original function and slightly modified the functional form for use with NonlinearModelFit.

nlm3D = NonlinearModelFit[z2, a + b Exp[-(x - mux)^2/(2 vx) - (y - muy)^2/(2 vy)],
   {{a, 25}, {b, Max[data]}, {mux, 100}, vx, {muy, 100}, vy}, {x, y}]
nlm3D["ParameterTable"]
estimates = nlm3D["BestFitParameters"]

with output

The following creates a visual assessment of the fit around the signal:

p1 = ListPointPlot3D[z2, 
   PlotRange -> {{80, 115}, {80, 115}, {0, 20000}}, 
   BoxRatios -> {1, 1, 0.5}, 
   AxesLabel -> Map[Style[#, Bold, Large] &, {"x", "y", "z"}], 
   PlotStyle -> PointSize[0.01]];
p2 = Plot3D[
   a + b Exp[-(x - mux)^2/(2 vx) - (y - muy)^2/(2 vy)] /. estimates, {x, 80, 115}, {y, 80, 115}, 
   PlotRange -> {{80, 115}, {80, 115}, {0, 20000}}, 
   PlotStyle -> Opacity[0.2]];
Show[{p1, p2}]

It sure seems like there is truncation for the higher values. Whether or not this is an adequate fit depends on your objective.

Your dataset has 202 rows and 101 columns of integers with about 47% of the values being zero. Are these frequency counts or measurements? I ask because the extremes in the marginal distributions have counts way too high to have come from a bivariate normal distribution. The marginal distributions (i.e., summing over the rows and columns) do have a "normal" shape in the middle but again if these are frequency counts, then those counts are way too high in the extremes to come from just a bivariate normal.

If you would describe the process that generates the data, then more constructive suggestions might be forthcoming. (Maybe what you have is a weighted mixture of a bivariate uniform distribution and a bivariate normal - if the values are frequency counts.)

Thnx for your nice code. For integrating the nlm3D, should I use Integrate command?

How would it be if I want to use FindFit for this matter?