OK, I "generalized" the Code a bit, so that it should work with other values for the experimentally fixed Ha and Da. See here and last part of the attachement.

Clear[t2]
t2[pK1_, pK2_, Ga_, Haa_, Daa_] := Evaluate[
  Eliminate[eqs /. {Ha -> Haa, Da -> Daa}, {H, G, D, HG, Ha}]]


Clear[f1]
f1[pK1_?NumericQ, pK2_?NumericQ, Ga_?NumericQ, Ha_?NumericQ, 
  Da_?NumericQ] := Module[{xx},
  xx = HD /. NSolve[t2[pK1, pK2, Ga, Ha, Da], HD];
  Part[Select[Flatten[xx], 0 < # < Ha &], 1]
  ]

In[14]:= FindFit[data, f1[a, b, x, 1, 1], {a, b}, x]

Out[14]= {a -> 2., b -> 3.}

But you should check it out with a couple of trials, meaning you should generate "experimental data" like your data above and check which answers are given by FindFit.

(And yes, years ago I worked together with Prof Hahn in a technical project and had some private contacts as well)

Attachments:

Well, your problem is somewhat mysterious because you don't say which are experimental inputs.

I assume Ha and Da are values which can be "fixed" experimentally. With this I do

eqs = {10^pK1 == HD/(H*D), 10^pK2 == HG/(H*G), Ha == HD + HG + H, Ga == HG + G, Da == HD + D};

Now I get rid of some parameters taking the experimentally fixed values into account:

In[6]:= t1 = Eliminate[eqs /. {Ha -> 1, Da -> 1}, {H, G, D, HG, Ha}]

Out[6]= -10^(2 pK1) + 10^pK1 HD + 3 10^(2 pK1) HD - 
  10^(pK1 + pK2) HD + 2^(1 + pK1 + pK2) 5^(pK1 + pK2) HD^2 - 
  10^pK1 HD^2 - 3 10^(2 pK1) HD^2 + 10^pK2 HD^2 + 10^(2 pK1) HD^3 - 
  10^(pK1 + pK2) HD^3 == 10^(pK1 + pK2) Ga (-1 + HD) HD

Then with t1 I define a function to get HD for a given Ga "knowing" that 0 < HD < 1 (because Ha == HD + HG + H and Ha == 1 )

Clear[ff]    
    ff[pK1_?NumericQ, pK2_?NumericQ, Ga_?NumericQ] := Module[{xx},
      xx = NSolve[
        -10^(2 pK1) + 10^pK1 HD + 3 10^(2 pK1) HD - 10^(pK1 + pK2) HD + 
          2^(1 + pK1 + pK2) 5^(pK1 + pK2) HD^2 - 10^pK1 HD^2 - 
          3 10^(2 pK1) HD^2 + 10^pK2 HD^2 + 10^(2 pK1) HD^3 - 
          10^(pK1 + pK2) HD^3 == 10^(pK1 + pK2) Ga (-1 + HD) HD, HD];
      Part[Select[Flatten[HD /. xx], 0 < # < 1 &], 1]
      ]

Test :
In[24]:= ff[2, 3, .6]

Out[24]= 0.455658

Then with the values of yours given above

data = {{1.*10^-8, 0.904875}, {0.05, 0.878427}, {0.1, 
    0.846593}, {0.15, 0.81062}, {0.2, 0.7719}, {0.25, 0.731591}, {0.3,
     0.690551}, {0.35, 0.649405}, {0.4, 0.608618}, {0.45, 
    0.568561}, {0.5, 0.529541}, {0.55, 0.491828}, {0.6, 
    0.455658}, {0.65, 0.42124}, {0.7, 0.38875}, {0.75, 
    0.358325}, {0.8, 0.330057}, {0.85, 0.30399}, {0.9, 
    0.280116}, {0.95, 0.25838}};


In[26]:= FindFit[data, ff[a, b, x], {a, b}, x]

Out[26]= {a -> 2., b -> 3.}

This should be checked with other combinations of parameters.

If possible (and you know him) give my greetings to Prof. Horst Hahn

Attachments:

Your code does not evaluate, because you do not provide the starting values of Ha etc. I am not sure what you want, but perhaps you could try a Manipulate:

hdPts[Ha_, Da_, Gmax_, pK1_, pK2_, displaypoints_] := 
  Module[{HD, HG, H, G, D},
   HDpoints = {}; HGpoints = {}; Hpoints = {}; Gpoints = {}; 
   Dpoints = {}; 
   Quiet[Flatten[
     Flatten[Delete[
       Do[{solut = 
          N[NSolve[{10^pK1 == HD/(H*D), 10^pK2 == HG/(H*G), 
             Ha == HD + HG + H, Ga == HG + G, Da == HD + D}, {HD, HG, 
             H, G, D}], 30], HDfirst = HD /. Part[solut, 1, 4], 
         HDsecond = HD /. Part[solut, 2, 4], 
         HDthird = HD /. Part[solut, 3, 4], 
         HGfirst = HG /. Part[solut, 1, 5], 
         HGsecond = HG /. Part[solut, 2, 5], 
         HGthird = HG /. Part[solut, 3, 5], 
         Hfirst = H /. Part[solut, 1, 1], 
         Hsecond = H /. Part[solut, 2, 1], 
         Hthird = H /. Part[solut, 3, 1], 
         Gfirst = G /. Part[solut, 1, 2], 
         Gsecond = G /. Part[solut, 2, 2], 
         Gthird = G /. Part[solut, 3, 2], 
         Dfirst = D /. Part[solut, 1, 3], 
         Dsecond = D /. Part[solut, 2, 3], 
         Dthird = D /. Part[solut, 3, 3],
         If[
          0 <= HDfirst <= Da && 0 <= HGfirst <= Ga && 
           0 <= Dfirst <= Da && 0 <= Gfirst <= Ga && 
           0 <= Hfirst <= Ha && HGfirst <= Ha && HDfirst <= Ha && 
           Hfirst + HGfirst + HDfirst <= Ha && 
           Gfirst + HGfirst <= Ga && 
           Dfirst + HDfirst <= Da, {AppendTo[
            HDpoints, {Ga, HDfirst/Ha}], 
           AppendTo[HGpoints, {Ga, HGfirst/Ha}], 
           AppendTo[Hpoints, {Ga, Hfirst/Ha}], 
           AppendTo[Dpoints, {Ga, Dfirst/Da}], 
           AppendTo[Gpoints, {Ga, Gfirst/Ga}]}],
         If[
          0 <= HDsecond <= Da && 0 <= HGsecond <= Ga && 
           0 <= Dsecond <= Da && 0 <= Gsecond <= Ga && 
           0 <= Hsecond <= Ha && HGsecond <= Ha && HDsecond <= Ha && 
           Hsecond + HGsecond + HDsecond <= Ha && 
           Gsecond + HGsecond <= Ga && 
           Dsecond + HDsecond <= Da, {AppendTo[
            HDpoints, {Ga, HDsecond/Ha}], 
           AppendTo[HGpoints, {Ga, HGsecond/Ha}], 
           AppendTo[Hpoints, {Ga, Hsecond/Ha}], 
           AppendTo[Dpoints, {Ga, Dsecond/Da}], 
           AppendTo[Gpoints, {Ga, Gsecond/Ga}]}],
         If[
          0 <= HDthird <= Da && 0 <= HGthird <= Ga && 
           0 <= Dthird <= Da && 0 <= Gthird <= Ga && 
           0 <= Hthird <= Ha && HGthird <= Ha && HDthird <= Ha && 
           Hthird + HGthird + HDthird <= Ha && 
           Gthird + HGthird <= Ga && 
           Dthird + HDthird <= Da, {AppendTo[
            HDpoints, {Ga, HDthird/Ha}], 
           AppendTo[HGpoints, {Ga, HGthird/Ha}], 
           AppendTo[Hpoints, {Ga, Hthird/Ha}], 
           AppendTo[Dpoints, {Ga, Dthird/Da}], 
           AppendTo[Gpoints, {Ga, Gthird/Ga}]}]}, {Ga, 0.00000001, 
         Gmax, Gmax/displaypoints}], 1], 1], 1]]];
Manipulate[
 ListLinePlot[hdPts[Ha, Da, Gmax, pK1, pK2, displaypoints], 
  PlotRange -> {0, 1}],
 {Ha, 0, 1},
 {Da, 0, 1},
 {Gmax, 0, 1},
 {pK1, 0, 1},
 {pK2, 0, 1},
 {{displaypoints, 100}, 1, 200}]

As it is, it does not work, perhaps because of bad initial values. It may not be a good idea to use capital letters as variables, specially one like D, which has its own meaning.