Thanks for the reply, so first I guess multiplying the covariance functions might also be possible The Kronecker Structured Covariance Matrix was probably the hint I needed. So in my case the gp.MarginalKron should hopefully work. Thank you very much!
Regdarding the stacking do I understand it rightly that you mean .eg. for y[0] = ( [y_dut1(t_0),…,y_dut1(t_n)] , t in T_1)? Then my question would be why the time steps are part of observations?
Maybe to clarify the aim here. I try to model a degrdation process under accelerated lifetime testing. At each stress condition S_i a set of J slightly different timesieres were recorded, let’s call them y_j = f(t,S_i) = exp(alpha_i * t). Here alpha_i is the decay rate. For each S_i alpha can be assumed as N(µ_i, sigma_i) (intersample variations) and with increased stress S_i-1 < S_i also follows µ_i-1, sigma_i-1 < µ_i, sigma_i (hingher stress means faster decay).
In this regard I have two follow up questions. So in the final application there will be about 5000 observations with a test structure of 4 different condition types (t, S1, S2, S3) yielding a y shape of (5000 x 1) and X with (5000x4). This would result in a covariance matrix in the dimensionality region of somewhat (1e15 x 1e15). Would it be a suitable approach (however that my work) to first model an individual GP for each f_j(t) at a given set of (S1,S2,S3) and finally model a GP of those individual GPs. This would drastically reduce the dimensionality of the covariance matrix and thus spead up the process. I hope you get my idea 
Thanks already for the input!