I have another problem with the aforementioned fit. I would like to perform the fit considering also as an independent variable "t". The latter is described by two similar equations, the first (T1) one for 0 <= t <= tau and the second (T2) for t > tau.

T2[x_, y_, t_, b_, l_, d_, tau_, m_, Intens_] /; 
  NumberQ[x] && NumberQ[y] && NumberQ[t] && NumberQ[b] && NumberQ[l] &&
    NumberQ[d] && NumberQ[tau] && NumberQ[m] && NumberQ[Intens] := 
 Intens/(2*\[Pi]*K0)*Sqrt[1 + m^2]*
  NIntegrate[
   Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\[Alpha] - 
            d))^2]/Sqrt[
       Dif*4*t]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\
\[Alpha] - d))^2]) - 
    Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\[Alpha] - 
            d))^2]/Sqrt[
       Dif*4*(t - 
          tau)]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\
\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha], 
    l*(1/(1 + m^2))}]

In producing noisy data I already get the first problems:

Block[{x = 0, y = 0, b = 0.001, l = 0.001, d = 0.0001, 
   m = -Tan[45*Pi/180], Intens = 25000}, 
  datat1 = Table[{0, 0, t, 
     T1[x, y, t, b, l, d, m, Intens] + 
      Random[NormalDistribution[0, 0.1]]}, {t, 0, 1, 0.01}]];

NIntegrate::inumr: The integrand (Sqrt[2] Erfc[ComplexInfinity])/Sqrt[\[Alpha]^2+(0.0001 +\[Alpha])^2+\[Beta]^2] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,0.0005},{0.,0.0005}}.
Power::infy: Infinite expression 1/0. encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.

Then I thought about doing the same thing for T2 and then using Join for the 4 data sets. But I think the problem starts already in the creation of noisy data.

I also tried with:

If[0 <= t <= 1, T1[0, 0, t, b, l, d, m, Intens], 
 T2[0, 0, t, b, l, d, 1, m, Intens]]

but nothing happens.

I wish I could fit the 4 datasets (for x, y and t) together with the 2 equations, always if this is possible.

You are fitting multivariate data (function dependent on x and y) but your data only has one coordinate logged (either x or y). so either add the respective x and y values to your data or just generate the full data range for x and y.

Both are equally valid in this case.

Dif = 0.00000013; K0 = 0.5; z = 0; linf\[Alpha] = 0.0;
T1[x_, y_, t_, b_, l_, d_, m_, Intens_] /; 
  NumberQ[x] && NumberQ[y] && NumberQ[b] && NumberQ[l] && NumberQ[d] &&
    NumberQ[m] && NumberQ[Intens] := 
 Intens/(2*\[Pi]*K0)*
  NIntegrate[
   Sqrt[1 + m^2]*
    Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\[Alpha] - 
              d))^2]/Sqrt[
        Dif*4*t]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\
\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha], 
    l*(1/(1 + m^2))}]

t = 1;
b = 0.001;
l = 0.001;
d = 0.0001;
m = -Tan[45*Pi/180];
Intens = 25000;

step = 0.00006;

datax = Table[{x, 0, 
    T1[x, 0, t, b, l, d, m, Intens] + 
     Random[NormalDistribution[0, 0.1]]}, {x, -0.002, 0.002, step}];
datay = Table[{0, y, 
    T1[0, y, t, b, l, d, m, Intens] + 
     Random[NormalDistribution[0, 0.1]]}, {y, -0.002, 0.002, step}];
alldata = Join[datax, datay];

dataxy = Flatten[
   Table[{x, y, 
     T1[x, y, t, b, l, d, m, Intens] + 
      Random[NormalDistribution[0, 0.1]]}, {x, -0.002, 0.002, 
     4 step}, {y, -0.002, 0.002, 4 step}], 1];

Show[
 Plot[{T1[x, 0, t, b, l, d, m, Intens], 
   T1[0, x, t, b, l, d, m, Intens]}, {x, -0.002, 0.002}, 
  PlotStyle -> Red],
 ListPlot[{datax[[All, {1, 3}]], datay[[All, {2, 3}]]}, 
  PlotStyle -> {Black, Gray}]
 ]

Show[
 Plot3D[
  T1[x, y, t, b, l, d, m, Intens], {x, -0.002, 0.002}, {y, -0.002, 
   0.002},
  PlotRange -> Full, Mesh -> False, PlotStyle -> Gray, 
  Lighting -> "Neutral"],
 ListPointPlot3D[dataxy, PlotRange -> Full, PlotStyle -> Black],
 Graphics3D[{Thick, Red, Line[datax], Line[datay]}]
 ]

(*only lines*)
fitL = FindFit[alldata, 
  T1[x, y, 1, bf, lf, d, m, Intens], {{bf, 0.005}, {lf, 0.005}(*,d,m,
   Intens*)}, {x, y}, Method -> "LevenbergMarquardt"]

(*full surface*)
fitS = FindFit[dataxy, 
  T1[x, y, 1, bf, lf, d, m, Intens], {{bf, 0.005}, {lf, 0.005}(*d,m,
   Intens*)}, {x, y}, Method -> "LevenbergMarquardt"]

That is interesting.

I had a look at your data ( datax and datay ) and tried this:

ff[x_, a_, b_, c_] := a  Exp[-b (x - c)^2]

and

sola =  FindFit[datax, ff[x, a, b, c], {{a, 5.85}, {b, 5.1 10^6}, {c, 0.00015}}, x]

(* {a -> 5.54423, b -> 4.52534*10^6, c -> 0.000158701} *)

Then

Plot[Evaluate[ff[x, a, b, c] /. sola], {x, -0.002, 0.002},
 Epilog -> Point /@ datax, PlotRange -> {0, 7}]

Similarly

solb =  FindFit[datay,   ff[x, a, b, c], {{a, 5.05}, {b, 2.56 10^6}, {c, -0.00001}}, x]
(* {a -> 5.54337, b -> 2.59886*10^6, c -> 1.4393*10^-7} *)
Plot[Evaluate[ff[x, a, b, c] /. solb], {x, -0.002, 0.002},
 Epilog -> Point /@ datay, PlotRange -> {0, 7}]

These are the results from FindFit, and I believe it. Anyhow, one could get the impression (optically) that there is somethig better for datax, but note, this is an impression. Play around with a, b and c:

Manipulate[
 Plot[ff[x, a, b, c], {x, -0.002, 0.002}, PlotRange -> {0, 7},
  Epilog -> Point /@ datax
  ],
 {a, 5, 7}, {b, .5 10^6., 10^7}, {c, -.002, .002}]

If I reformulate with Integrate instead of NIntegrate it does not give an error:

FindFit[alldata,
 0.3183098861837907` Intens Integrate[
   Sqrt[1 + m^2]*
    Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\[Alpha] - 
              d))^2]/Sqrt[
        Dif*4*1]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\
\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha], 
    l*(1/(1 + m^2))}],
 {b, l, d, m, Intens}, {x, y}, Method -> "LevenbergMarquardt"]

but it takes such a long time that I did not wait. The numerical evaluation of the double integral for so many different values of the parameters may be very onerous.