Hi,

I am new to ASPECT and have a basic confusion while developing a depth-dependent material model.

I have a 4000 × 660 km 2-D box model, where I am interested in implementing different viscosity layers. For simplicity, I am keeping 2 layers: 0-100 km and 100-660 km. I want the top layer to have a viscosity of 1e21 Pa·s and the bottom layer to have a viscosity of 1e20 Pa·s. I am preparing the PRM file as follows:

```
subsection Material model
set Model name = depth dependent
subsection Depth dependent model
set Base model = simple
set Reference viscosity = 1e21
set Depth dependence method = List
set Depth list = 100.0000000e+00, 6.7000000e+05
set Viscosity list = 1.0000000e+21, 1.0000000e+20
end
subsection Simple model
set Thermal conductivity = 1e-6
set Thermal expansion coefficient = 3e-5
set Reference density = 3300
set Reference temperature = 1350
set Density differential for compositional field 1 = 1
set Composition viscosity prefactor = 1
end
end
```

The resulting viscosity becomes uniform at approximately 5e23 Pa·s, and the model runs for only a few time steps (for this particular setup, it took only 6 steps) to reach 3.5 billion years.

I have tried different values for the reference viscosity and viscosity list but have never obtained the correct result.

According to the documentation and source code: [ \eta = \eta(z) * \eta(T, X, P, \ldots) / \eta(r) ]

But the result in my case looks very different. So, I am not clear on how the net viscosity is calculated using this method. Any explanation would be helpful.