Hi Francis,

Could you please give me some instructions on how to fix them [`μs`

, `qDy`

] to accept the heterogenous background model?

For the shear viscosity, you could to following:

- comment out line 114;
- define a “host” array for
`μs`

: `Mus = zeros(nx ,ny )`

after line 145;
- initialise
`Mus`

with your newly implemented spatially variable background porosity array `Phi0`

, or with another function of `Phi`

, e.g. `Mus .= ηC0*Phi0/η2μs`

;
- copy the host array to the device array
`Mus = Data.Array(Mus)`

.

Then one will have to modify all the places where the scalar (`Data.Number`

) `μs`

was used, replacing it by the `Data.Array`

array `Mus`

. E.g. kernel definitions will have to take `Data.Array::Mus`

as argument now.

Since the shear viscosity is now an array, you will have to change the definition of the pseudo-time steps. You can to comment out lines 120-121, initialise 3 array:

```
dτVx = @zeros(nx-1,ny-2)
dτVy = @zeros(nx-2,ny-1)
dτPt = @zeros(nx ,ny )
```

and the following kernel to compute them:

```
@parallel function compute_dτ!(...)
@all(dτVx) = min_dxy2/@av_xi(Mus)/(1.0+β_n)/4.1/Vsc
@all(dτVy) = min_dxy2/@av_yi(Mus)/(1.0+β_n)/4.1/Vsc
@all(dτPt) = 4.1*@all(Mus)*(1.0+β_n)/max(nx,ny)/Ptsc
return
end
```

taking care to pass the correct arguments to the new kernel (here `...`

). See the 2D Stokes miniapp for reference.

The y-direction Darcy flux `qDy`

calculation you point to is the vertical boundary condition definition of the inflow and outflow of the model. The definition is chosen to satisfy an effective pressure `Pt-Pf=0`

at the boundaries to avoid compaction/decompaction at the boundaries. If you want to have non-scalar values there, you could perform vectorised calculation on host arrays or create a kernel to precompute a vector for the lower and upper boundary. This may actually not be really needed.

Could you also give me some brief instruction on adapting time stepping […] how should I change these time stepping

There is not much to do there. You can play with `dt_red`

, the controls on the physical time step reduction, to see if you better resolve the waves depending on the configuration.

Also, if we change the porosity to be heterogenous, permeability will also be heterogenous. How do we calculate the compaction length then?

Porosity is heterogenous. What you change is the background porosity, the “reference value” of the variable porosity. Permeability is and will further be heterogenous as it depends on `Phi`

and not on on `ϕ0`

. To compute the compaction length, you can choose which reference permeability you want to consider. I would take the far field, background (or reference) one.

I hope this helps. If you need further support, feel free to start a dev git repo; I could then provide more direct support for your particular needs and project.

Ludovic