I think you are making a silly mistake by considering solution as a List when you have defined it as a InterpolatingFunction object. Here is little modification which may help
Clear[Eq0, FilmEqn, h, Bo, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\) +
Div[-h^3 Bo Grad[h] +
h^3 Grad[Laplacian[h]] + (\[Delta] h^3)/(Bi h + K1)^3 Grad[h] +
m (h/(K1 + Bi h))^2 Grad[h]] + \[Epsilon]/(
Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x];
SetCoordinates[Cartesian[x, y, z]];
FilmEqn[Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_, time_] :=
Eq0[solution[x, y, time], {Bo, \[Epsilon], K1, \[Delta], Bi, m, r}]
expr = FilmEqn[0, 10^-6, 1, 10^-3, 1, 0.05, 0, time];
fun[a_, b_, t_] := Evaluate[expr /. x -> a /. y -> b /. time -> t];
Plot3D[fun[x, y, 100000], {x, 0, L}, {y, 0, L}, PlotPoints -> 40,Mesh ->None,
ColorFunction -> Hue, PlotRange -> All]
