mbakker7 / timml

An analytic element model for steady multi-layer flow
MIT License
38 stars 21 forks source link

Inconsistent results when Polygon Inhomogenity (Semi-Confined Top) in combination with a Leaky Line Doublet #63

Closed ArtemisRo closed 1 year ago

ArtemisRo commented 1 year ago

I am working on a groundwater model in Qgis-Tim which includes a Building Pit. However, the Building Pit did not work as expected in Qgis-Tim, so I replaced it with a Polygon Inhomogenity (Semi-Confined Top) in combination with a Leaky Line Doublet.

I have encountered an issue with the use of Leaky Line Doublet and Polygon Semi-Confined Top elements together. Specifically, when the size of the Polygon Inhomogenity element is exactly the same as the Leaky Line Doublet, the calculated head in layer 1 becomes extremely high. I have attached a figure with the contour lines to illustrate the issue. Furthermore, I added the script below:

As a workaround, I made the Polygon Inhomogenity slightly larger than the Leaky Line Doublet and then the results are reasonable.

`import numpy as np import timml

model = timml.ModelMaq( kaq=[5.0, 19.0, 19.0, 19.0, 19.0, 40.0], z=[ -3.5, -3.5, -5.5, -7.8, -16.0, -16.1, -32.0, -32.1, -64.0, -65.0, -71.0, -72.0, -185.0, ], c=[1000.0, 1200.0, 6.0, 12.0, 25.0, 45.0], topboundary="semi", npor=[None, None, None, None, None, None, None, None, None, None, None, None], hstar=-4.3, )

damwand1_0 = timml.LeakyLineDoubletString( xy=[ [125346.50031011457, 479033.093129771], [125381.72290076336, 479053.37175572524], [125389.83435114504, 479038.0706106871], [125353.51717557252, 479018.8980916031], [125346.50031011457, 479033.093129771], ], res=100.0, layers=0, order=4, label="dw0", model=model, ) damwand1_1 = timml.LeakyLineDoubletString( xy=[ [125346.50031011457, 479033.093129771], [125381.72290076336, 479053.37175572524], [125389.83435114504, 479038.0706106871], [125353.51717557252, 479018.8980916031], [125346.50031011457, 479033.093129771], ], res=100.0, layers=1, order=4, label="dw1", model=model, )

injectie3_0 = timml.PolygonInhomMaq( kaq=[5.0, 19.0, 19.0, 19.0, 19.0, 40.0], z=[ -3.5, -3.5, -5.5, -7.8, -16.0, -16.1, -32.0, -32.1, -64.0, -65.0, -71.0, -72.0, -185.0, ], c=[50.0, 1200.0, 6.0, 12.0, 25.0, 45.0], topboundary="semi", npor=[None, None, None, None, None, None, None, None, None, None, None, None], hstar=-4.3, xy=[ [125346.50031011457, 479033.093129771], [125381.72290076336, 479053.37175572524], [125389.83435114504, 479038.0706106871], [125353.51717557252, 479018.8980916031], [125346.50031011457, 479033.093129771], ], order=4, ndeg=6, model=model, ) Headwell1_0 = timml.HeadWell( xw=125370.84029240627, yw=479032.82079650735, hw=-8.3, rw=0.3, res=1.0, layers=1, label="w1", model=model, ) model.solve() head = model.headgrid( xg=np.arange(125330.0, 125427.0, 2), yg=np.arange(479071.0, 478994.0, -2) )`

Model results

MattBrst commented 1 year ago

I would suggest to use the min/max head as diagnostic information (for example in QGIS-TIM). In the example above, a head of +200 m is calculated, which is not very realistic.

dbrakenhoff commented 1 year ago

Perhaps we can add some check for control points instead of checking heads? The heads check is challenging to do in a robust way I think.

As stated by @rt84ro, the above problem is fixed when nudging the LeakyLineDoublet by a small distance, e.g. 1e-5:

eps = 1e-5
xy = np.array(
    [
        [125346.50031011457 - eps, 479033.093129771 + eps],
        [125381.72290076336 + eps, 479053.37175572524 + eps],
        [125389.83435114504 + eps, 479038.0706106871 - eps],
        [125353.51717557252 - eps, 479018.8980916031 - eps],
        [125346.50031011457 - eps, 479033.093129771 + eps],
    ]
)

Overlapping control points are not necessarily a problem though. Inhomogeneities add overlapping points for each side (inhom.Nsides * (inhom.order+1)), and overlapping points from different elements in different layers are also fine I believe. So maybe that's also difficult to implement in a robust way.

mbakker7 commented 1 year ago

The building pit element was developed especially to combine an inhomogeneity and an impermeable wall. This works well in TimML. The implementation in QGIS-TIM apparently included a bug but that will be fixed soon. TimML will be modified to make the boundary of the building pit leaky rather than impermeable. We should include in the docs that these elements can not be combined.