Issues with using certain coordiante values with particle tracking

[waterrocks247]

I'm trying to use the particle tracking code in ttim to look at ASR recovery volumes and I'm getting some weird results depending on the coordinates I use for the particles and wells. If I use coordinates in feet or meters from a CRS as inputs for well locations and thus particle locations, the particles all track shoot off into one direction. However using a 0 based coordinate system where the well and particles begin from 0,0 the particles will radiate symetrically from the well. Also the particles will begin to radiate more symetrically from the well the greater distance away I place them from the well. I thought this behavior was quite strange so I provided code below to reporduce what I am seeing.

import numpy as np
import matplotlib.pyplot as plt
import timml as tm
import ttim as tt
# parameters
Q = 100 # discharge of well, m^3/d
k = 10 # hydraulic conductivity, m/d
H = 10 # thickness of aquifer, m
Ss = 1e-4 # specific storage, m^(-1)
npor = 0.3 # porosity, -
xw = 0 # x-location of well
yw = 0 # y-location of well
rw = 0.3 # radius of well, m
tmin = 0.001 # first time of simulation after change in bc, d
ml = tt.ModelMaq(kaq=[k], z=[H, 0], Saq=[Ss], poraq=npor,
tmin=tmin, tmax=1000, M=10)
crs = True
if crs:
    w = tt.Well(ml, xw=5688794 , yw=19649311, tsandQ=[(0, Q), (10, -Q)], rw=rw)
else:
    w = tt.Well(ml, xw=0 , yw=0, tsandQ=[(0, Q), (10, -Q)], rw=rw)
ml.solve()
alpha = np.linspace(0, 2 * np.pi, 25)
# y: 19649311.549242355
# X: 5688794.751554074 
if crs:
    xt =5688794+ rw * np.cos(alpha)
    yt =19649311 + rw * np.sin(alpha)
else:
    xt = rw * np.cos(alpha)
    yt = rw * np.sin(alpha)
xyztend = []
xyzts =[]
plt.subplot()
for i in range(len(xt)):
    trace = tt.timtrace(ml, xstart=xt[i], ystart=yt[i], zstart=0.5 * H, tstartend=[0, 10,20],tstartoffset=tmin, tstep=1, hstepmax=2,silent=True)
    xyzt = trace['xyzt']

    print('xyzt:')
    print(trace['xyzt'])
    # if trace['status'][0] < 0: # likely undesirable
    #     print(trace['message'])
    plt.plot(xyzt[:, 0], xyzt[:, 1])
    xyztend.append(xyzt[-1])
    xyzts.append(xyzt)
plt.grid()
xyztend = np.array(xyztend)
plt.plot(xyztend[:, 0], xyztend[:, 1], 'k--', label='not recovered')
plt.scatter(xyztend[:,0], xyztend[:,1], marker='o', label = 'End Position')
plt.legend()
plt.xlabel('x (m)')
plt.ylabel('y (m)');

plt.figure(figsize=(10, 6))
for data in xyzts:
    # Separate columns for x, y, z, and time
    x = data[:, 0]
    y = data[:, 1]
    z = data[:, 2]  # This isn't used in this plot but could be useful for 3D plotting
    time = data[:, 3]
    
    # Create scatter plot with x, y, and color coded by time
    
    sc = plt.scatter(x, y, c=time, cmap='viridis', edgecolor='k')
plt.scatter(0,0, label='Well')
plt.colorbar(sc, label='Time')  # Add colorbar to show the time scale
# Labels and title
plt.xlabel('X Position')
plt.ylabel('Y Position')
plt.title('Position of Points with Color Coded by Time')
plt.grid(True)
plt.show()

[mbakker7]

Thanks for your question, @waterrocks247
I am not sure I understand the problem yet that you are trying to solve.
You specify tsandQ=[(0, Q), (10, -Q)].
This means that from t=0 till t=10 you are extracting water at a rate Q=100 m^3/d.
Then for t>10, you are injecting water at a rate Q=100 m^3/d.
Is that what you intended to do, first extracting and the injecting?

[waterrocks247]

I did actually intend to inject then extract, but the problem I am concerned with is what the particle tracks do given the same inputs but only change the location of the well and the particles. If I use the crs coordiantes the particles take a completely different course than not using the crs coordinates when nothing else has changed. So if you swap out the crs=True or False the particles have completely different behavior. This doesn't make a lot of sense to me because the x,y position of the well and particles shouldn't impact the results if the change was relative to 0.

[mbakker7]

You found an undesired feature of the pathline tracing code. The current implementation moves the specified starting location a very small amount to avoid getting stuck in a stagnation point in symmetric cases (like yours). The way this is programmed, this "very small amount" turns out to be significant in case the coordinates are very large. This should be fixed. Can you file an Issue?
I guess the easiest fix is probably be to make it an option to move the starting location a small amount. That way the user is in control.
Thanks for reporting this undesired behavior!

[waterrocks247]

Yes I can file an issue. It sounds like maybe a quick work around for me, in the case of mulitple injection wells is to make the coordinates of one well the origin at 0,0 and use this to transform the coordinates of other wells and that well so the numbers are small then I when I export the tracks to gis I can apply the transform to all the x y positions.

[mbakker7]

Yes, the work-around for now is: pick an origin somewhere in your model domain, or even right at your well. Translate all your coordinates and do your traces in this translated coordinate system.