I am trying to simulate a dataset that matches empirical data that I have, but am unsure how to estimate the errors in the original data. The empirical data includes heteroscedasticity, but I am not interested in transforming it away, but rather using a linear model with an error term to reproduce simulations of the empirical data.

For example, let's say I have some an empirical dataset and a model:

n=rep(1:100,2)
a=0
b = 1
sigma2 = n^1.3
eps = rnorm(n,mean=0,sd=sqrt(sigma2))
y=a+b*n + eps
mod <- lm(y ~ n)

using plot(n,y) we get the following.

However, if I try to simulate the data, simulate(mod), the heteroscedasticity is removed and not captured by the model.

I can use a generalized least squares model

VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)

that provides a better model fit based on AIC, but I don't know how to simulate data using the output.

My question is, how do I create a model that will allow me to simulate data to match the original, empirical data (n and y above). Specifically, I need a way to estimate sigma2, the error, using either using a model?