The biomass content and carbon captured by forest plantations is of interest, for example in the context of climate change and carbon budgets.The main objective of our study was to develop functions to estimate the total biomass and its components (stem, branches, bark and leaves) of Pinus radiata D. Don trees in Chile. The methodology proposed for the model fitting uses the maximum likelihood method in a multivariate equation system fitting simultaneously. The fit strategy incorporates additivity restrictions in the estimation functions and in the variance functions to incorporate the heteroskedasticity of biomass, and three structures of the variance-covariance matrix were evaluated to assess the dependence of the different components of tree biomass. Non-linear biomass functions that used the variable (DH)-H-2 performed best according to several indicators of goodness-of-fit (log-likelihood, Akaike Information Criterion and Bayesian Information Criterion) and estimation precision (root mean square error (RMSE), Bias and EI). The simple structure of both biomass and variance estimation functions was beta(1)((DH)-H-2)beta(2), and in the modelling system for total tree biomass RMSE between 54.1-54.4 kg (28-36%) were obtained. The three variance-covariance matrix structures evaluated did not generate clear differences in relation to the RMSE, bias and Error Index indicators. The structure of the variance-covariance matrix that incorporated explicitly in the system equations allowed modelling of the relationship between biomass components.