multivariate mean? [edit : because my question was ambiguous, I decided to rewrite it entirely, with some simplification but a lot more details on the experimental design]
Four independent 10m*10m plots each received sewage sludge from the same water treatment facility (this is precisely where pseudoplication occurs).
After a year of waiting, one "representative" 1kg soil sample was sampled from each plot.
Ten 1g-sub-samples were sampled from each soil sample, then suspended in 10ml water and agitated for 2 days in order to "equilibrate".
Only then, 1 (arbitrary) unit of labeled molecule was added into each soil suspension.
We want to follow the decrease of labeled molecule quantity as time goes on.
After 1 min of agitation, remaining labeled molecule in water of one soil suspension was measured destructively, giving a value "v" between 0 and 1. The value of w = 1-v was recorded.
Same after 2, 3, 4, 5, 6, 7, 8, 9 and 10 minutes of agitation.
So far, we have recorded 10 w values for each soil sample.
plot of w with time of agitation :

Because v decreases with time, w is increasing with time and can be efficiently modeled by the following model equation :
      w ~ a * time^b

This model was fit onto each of the 4 sets of 10 w values, giving four sets of (a,b) parameters.
My problem is : How to calculate / estimate a set of global (a,b) parameters ?
I think calculating average of each parameter (=> a_mean and b_mean) is not right. Fitting a model on all pooled 40 w values is no more right.
NB : As a side note, if I had one, and only one, parameter of interest (for example, a maximum value, or a mean) by soil sample, one good option would be to average those four values into one global average.
Thanks !
 A: Your model assumption is
$$
w = a ~ t^b
$$
with some small random variations.
Take the logarithm of that equation:
$$
\log w = \log a + b ~ \log t.
$$
You now have a linear relationship between $\log w$ and $\log t$ with two parameters, $\log a$ and $b$. Determining the parameters is a standard linear regression problem, and if for simplicity we assume normally distributed errors with constant variance, the parameters can be estimated by ordinary least squares.

I don't know which software you use, but in Matlab you would do the following: Define the time points:
t = (1 : 10)';

Let's assume that your $w$-values are contained in a matrix W with 10 rows for the different time points and 4 columns for the different samples. Set up a design matrix where the first column consists of repetitions of $\log t$ and the second is constant 1
X = [repmat(log(t), size(W, 2), 1), ones(numel(W), 1)];

and construct a vector of logarithmized data
y = log(W(:));

Now you can perform the linear regression using Matlab's backslash operator:
parest = X \ y;

The estimates of the parameters are
b = parest(1);
a = exp(parest(2));

A: Thanks to some simple R simulations, I came up with a satisfying answer to my own question.
The question asks what is the right way of estimating the true / global  set of parameters of a given model, given the model has been fitted to N different datasets, yielding N  vectors of parameters.
The answer is : calculate vector average from these N vectors.
How I came up with this answer that contradicts what I though at first is quite well demonstrated by a simulation (= backing up my answer with personal experience) :
Let's imagine we observe two variables X (independent variable), and Y (dependent variable). Assuming Y is linearly dependent on X, we can construct the following linear model :
 Y = a * X + b

Assume we know by construction that a_TRUE = 2 and b_TRUE = 10.
In the case of a and b being orthogonal to each other, one can draw a sample of N "a" values, and another sample of N "b" values.
 a_sample <- rnorm(10,a_TRUE,1)
 b_sample <- rnorm(10,b_TRUE,2)

here is one possible plot of the true equation (black solid) and all 10 individual equations (dashed red) :

By construction, a_TRUE can be estimated by averaging all individual vectors.
Here (sorry, I forgot to set seed so it is not reproducible), b_estimate = 9.7 +- 0.7 (95% CI), and a = 1.5 +- 0.8
By analogy, in the nonlinear case, it seems OK to average over all 4 (a,b) vectors.
However, with only four values to average by parameter, it is likely that 95% confident intervals will be huge.
I don't know, thought, if this is still OK when the parameters are correlated.
Comments are welcome.
