Hi,

that's a numerical issue, I believe. You have qualities with vastly different orders of magnitudes. That is numerically very challenging. Try this:

(*Parameters*)
k2 = (1.8*10^-14);
k4 = (6*10^-34);
k6 = (2.9*10^-11);
k7 = (2.2*10^-10);
k8 = (26.3*10^-12);
k9 = (8.9*10^-12);
k10 = (1.9*10^-15);
k13 = (9.37*10^-12);
k14 = (8.6*10^-12);
k15 = (2*10^-15);
k16 = (4.4*10^-31);
k17 = (7*10^-11);
k18 = (5.6*10^-12);

(*ICs*)

M = 2.5*10^19;
O2 = 0.21*(2.5*10^19);
H2O = 3.81*10^17;
NO = (10*10^-9)*(2.5*10^19);
NO2 = (1*10^-9)*(2.5*10^19);
RH = (100*10^-9)*(2.5*10^19);

(*Equations*)

solution = NDSolve[{\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(a[
      t]\)\) == -a[t]*((8.9*10^-3)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}]) + 
     k2*b[t]*d[t] + k9*j[t]*b[t] + k14*m[t]*b[t] - 
     k16*h[t]*a[t]*M, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(b[
      t]\)\) == ((8.9*10^-3)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*a[t] - 
     k2*b[t]*d[t] - k9*j[t]*b[t] - k14*m[t]*b[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(c[
      t]\)\) == ((8.9*10^-3)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*
      a[t] + ((3.65*10^-4)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*d[t] - 
     k4*c[t]*O2*M + k6*fo[t]*M, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(d[t]\)\) == -k2*b[t]*
      d[t] - ((3.65*10^-4)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*d[t] + 
     k4*c[t]*O2*
      M - ((3.65*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*d[t] - 
     k15*m[t]*d[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(fo[
      t]\)\) == ((3.65*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*d[t] - 
     k6*fo[t]*M - k7*fo[t]*H2O, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(h[t]\)\) == 
    2*k7*fo[t]*H2O - k8*i[t]*h[t] - k13*l[t]*h[t] + k14*m[t]*b[t] + 
     k15*m[t]*d[t] - k16*h[t]*a[t]*M, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(i[t]\)\) == -k8*i[t]*h[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(j[t]\)\) == 
    k8*i[t]*h[t] - k9*j[t]*b[t] - k18*k[t]*m[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(k[t]\)\) == 
    k9*j[t]*b[t] - k10*k[t]*O2, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(l[t]\)\) == 
    k10*k[t]*
      O2 - ((2.96*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*
      l[t] - ((4.25*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*l[t] - 
     k13*l[t]*h[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(m[t]\)\) == 
    k10*k[t]*O2 + 
     2*((2.96*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*l[t] + 
     k13*l[t]*h[t] - k14*m[t]*b[t] - k15*m[t]*d[t] - 2*k17*(m[t])^2 - 
     k18*j[t]*m[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(n[
      t]\)\) == ((2.96*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*
      l[t] + ((4.25*10^-5)*
        Piecewise[{{0.001, 
           t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
           21600 <= t <= 64800}, {0.001, 64800 < t}}])*l[t] + 
     k13*l[t]*h[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(o[
      t]\)\) == ((4.25*10^-5)*
       Piecewise[{{0.001, 
          t < 21600}, {Sin[((t - 21600)/43200)*Pi] + 0.001, 
          21600 <= t <= 64800}, {0.001, 64800 < t}}])*l[t], \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(p[t]\)\) == 
    k16*h[t]*a[t]*M, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(q[t]\)\) == 
    k17*(m[t])^2, \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(r[t]\)\) == k18*j[t]*m[t], 
   a[0] == NO2, b[0] == NO, c[0] == 0, d[0] == 0, fo[0] == 0, 
   h[0] == 0, i[0] == RH, j[0] == 0, k[0] == 0, l[0] == 0, m[0] == 0, 
   n[0] == 0, o[0] == 0, p[0] == 0, q[0] == 0, r[0] == 0}, {a, b, c, 
   d, fo, h, i, j, k, l, m, n, o, p, q, r}, {t, 0, 86400}, 
  WorkingPrecision -> 35, MaxSteps -> Infinity]

It will give an (expected) error message which you can ignore. It is a much more benign error message than the one that your original code gives.

If you plot this:

Plot[Evaluate[{a[t], b[t], c[t], d[t], fo[t], h[t], i[t], j[t], k[t], l[t], m[t], n[t], o[t], p[t], q[t], r[t]} /. solution], {t, 0, 86400}, PlotRange -> All]

everything stays nice and positive. With this precision and MaxSteps value, you can then also integrate for longer times, e.g. up to 10000000. The thing is that below 80000 the system has (practically) reached its fixed point and nothing interesting happens after that.

Cheers,

Marco

PS: You might want to use:

Manipulate[Plot[Evaluate[s /. solution], {t, 0, 86400}, PlotRange -> All], {s, {a[t], b[t], c[t], d[t], fo[t], h[t], i[t], j[t], k[t], l[t], m[t], n[t], o[t], p[t], q[t], r[t]}}]

To look at the individual functions. In particular you might want to look at $h(t)$ which drops to very low values at times of around 65000 and has a relatively large negative slope just before that. Small numerical errors might make this go negative and lead to all sorts of trouble.