Ground based loop TDEM amplitude anomalies

Hi all,

As a continuation of my efforts to simulate a 3D TDEM survey with a source loop and a grid of receivers, I have come up with the following code, which I’ll post and annotate.
The code runs without errors, but the dB/dt values in the final output are clearly impossible (~e123). Here is what I am doing.
-I have a “filament-like” feature, depicted by a highly conductive dipping tube in an otherwise high resistivity model.
-The source is a ground based (z=0) loop with a radius of ~28 meters.
-Receivers are a uniform grid also on the surface, and record at 31 logspaced points between 1e-7 and 1e-3 (though I have tried various other time ranges).
-I am using a ramp on ramp off style of waveform (I tried others also, always with anomalous amplitudes as a result).
Any help in resolving this would be much appreciated. Here’s the code in blocks:

I am using a few extra packages, namely plotly to visualize the survey.

from discretize import TensorMesh
from discretize import TreeMesh
from discretize.utils import mkvc, refine_tree_xyz
import matplotlib.pyplot as plt
from SimPEG.utils import plot2Ddata, surface2ind_topo
from SimPEG import maps
import SimPEG.electromagnetics.time_domain as tdem
from SimPEG.utils import meshTensor
import numpy as np
import numpy.matlib as matlib
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
import plotly.offline as pyo
import plotly.graph_objs as go
import plotly.express as px
import plotly.io as pio
pio.renderers.default='browser'

%matplotlib inline
try:
    from pymatsolver import Pardiso as Solver
except ImportError:
    from SimPEG import SolverLU as Solver

save_file = False

Set up the mesh. I’m using a tensormesh here, with padding. The plotting is commented out. The last line of this section involves the model mapping, which I find somewhat opaque to understand (it seems to have to do with flagging the active cells in a mesh, but in my case ind_active returns all True, since it is not an air-based TDEM study. Is this correct?)

nc = 40  # number of core mesh cells in x, y and z
ncz = 24
dh = 10  # base cell width in x, y and z
dhz = 5
npad = 5  # number of padding cells: these are cells of varying size OUTSIDE the code cell region (usually less dense)
exp = 1.25  # expansion rate of padding cells

hx = [(dh, nc)] #reads like this: dilation left of zero, core cells, dilation positive from zero
hy = [(dh, nc)]
hz = [(dhz, 4, -exp), (dhz, ncz)]
mesh = TensorMesh([hx, hy, hz], x0="CCN") # creates the mesh reference

# The bottom southwest corner
x0 = mesh.x0 #indicates the southwest deepest cell

# The total number of cells
nC = mesh.nC

# An (nC, 3) array containing the cell-center locations
cc = mesh.gridCC

# A boolean array specifying which cells lie on the boundary
bInd = mesh.cellBoundaryInd

# The cell volumes
v = mesh.vol

# Plot all cells volumes or plot cell volumes for a particular horizontal slice
# fig = plt.figure(figsize=(9, 4))
# ax1 = fig.add_subplot(121)
# ax2 = fig.add_subplot(122)

# mesh.plotImage(np.log10(v), grid=True, ax=ax1)
# ax1.set_title("All Cell Log-Volumes")

# cplot = mesh.plotSlice(np.log10(v), grid=True, ax=ax2, normal="Z", ind=2)
# cplot[0].set_clim(np.min(np.log10(v)), np.max(np.log10(v)))
# ax2.set_title("Cell Log-Volumes #2")

air_value = 1e-8
ind_active = mesh.gridCC[:, 2] < 0.0
model_map = maps.InjectActiveCells(mesh, ind_active, air_value)

Define the conductivity model. sig_silver is the conductivity of the dipping tube feature.

sig_silver     = 10000
sig_background = 2e-3

model = np.log(sig_background) * np.ones(mesh.nC) #background model

#create line in 3D
xo,yo,zo = (0,0,-20) #passes through this point
x1,y1,z1 = (400,400,-250)
Lx,Ly,Lz = [np.linspace(xo,x1,100),np.linspace(yo,y1,100),np.linspace(zo,z1,100)]

#function to calculate the distance of every cell to a point
def dist(x,msh):
    dis = np.zeros(len(msh))
    p1 = np.array([mesh.gridCC[:,0], mesh.gridCC[:,1], mesh.gridCC[:,2]])
    p2 = np.transpose(matlib.repmat(np.transpose(x),len(p1[0]),1))
    dist = (((p2[0,:]-p1[0,:])**2)+((p2[1,:]-p1[1,:])**2)+((p2[2,:]-p1[2,:])**2))**(1/2)
    return dist

#Note: mesh.gridCC[x1,xyz] returns one element (x) and its 3d x,y,z position in the volume
thk = 10 #vein thickness

for i in range(len(Lx)):
    ds = dist(np.array([Lx[i],Ly[i],Lz[i]]),mesh)
    model[ds < thk] = np.log(sig_silver)

Define the receiver setup. Simple uniform grid of surface receivers of the dB/dt type.

# Observation times for response (time channels)
n_times = 31
time_channels = np.logspace(-7, -3, n_times)

n_tx = 11
# Define receiver locations
xrx, yrx, zrx = np.meshgrid(
    np.linspace(-190, 190, n_tx), np.linspace(-190, 190, n_tx), [0]
)
receiver_locations = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)]
receivers_list = []

ntx = np.size(xrx)
# Each unique location defines a new receiver
for ii in range(ntx):
    # Here we define receivers that measure the h-field in A/m
    dbzdt_receiver = tdem.receivers.PointMagneticFluxTimeDerivative(
        receiver_locations[ii, :], time_channels, orientation = "z"
    )
    receivers_list.append(dbzdt_receiver)

Source waveform is a step on, step off waveform, with a stable on time plateau and a final off time of 0s.

waveform_times = np.linspace(-0.002, 0, 21)
waveform = tdem.sources.TrapezoidWaveform(
    ramp_on=np.r_[-0.007, -0.006], ramp_off=np.r_[-0.001, 0.0], offTime=0.0
)

waveform_value = [waveform.eval(t) for t in waveform_times]

For those who want visualize the survey, here it the plotly script I’m using.

fig =go.Figure(data=go.Isosurface(
    x= mesh.gridCC[:,0],
    y= mesh.gridCC[:,1],
    z= mesh.gridCC[:,2],
    value=model,
    isomin=-6,
    isomax=10,
))

fig.add_scatter3d(    
    x=receiver_locations[:,0],
    y=receiver_locations[:,1],
    z=receiver_locations[:,2],
    mode='markers',
)

fig.show()

Here is the source layout as a loop with a radius of 28m and a current of 2.5A centered at [0,0,0].

src_ramp_on = tdem.Src.CircularLoop(
    receivers_list,
    location=np.r_[0.0, 0.0, 0.0],
    orientation = "z",
    current = 2.5,
    radius=28,
    waveform=waveform,
)
src_list_ramp_on = [src_ramp_on]

survey_ramp_on = tdem.Survey(src_list_ramp_on)

survey = tdem.Survey(src_list_ramp_on)

Run the simulation.

time_steps = [(1e-5, 20), (1e-5, 10), (1e-4, 30), (1e-3, 5)]

#check total time
dt = meshTensor(time_steps)
tmax = np.hstack([[0], dt]).cumsum()

#Note: waveform on time was defined as -0.007 to 0
simulation = tdem.simulation.Simulation3DMagneticFluxDensity(
    mesh, survey=survey, sigmaMap=model_map,verbose=True, solver=Solver, t0=-0.007
)

# Set the time-stepping for the simulation
simulation.time_steps = time_steps

dpred = simulation.dpred(model)

Plot the first and last times of the recordings.

# Data were organized by location, then by time channel
dpred_plotting = np.reshape(dpred, (n_tx ** 2, n_times))

# Plot
fig = plt.figure(figsize=(10, 4))

# dB/dt at early time
v_max = np.max(np.abs(dpred_plotting[:, 0]))
ax11 = fig.add_axes([0.05, 0.05, 0.35, 0.9])
plot2Ddata(
    receiver_locations[:, 0:2],
    dpred_plotting[:, 0],
    ax=ax11,
    ncontour=30,
    clim=(-v_max, v_max),
    contourOpts={"cmap": "bwr"},
)
ax11.set_title("dBz/dt at 0.0000001 s")

ax12 = fig.add_axes([0.42, 0.05, 0.02, 0.9])
norm1 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)
cbar1 = mpl.colorbar.ColorbarBase(
    ax12, norm=norm1, orientation="vertical", cmap=mpl.cm.bwr
)
cbar1.set_label("$T/s$", rotation=270, labelpad=15, size=12)

# dB/dt at later times
v_max = np.max(np.abs(dpred_plotting[:, -1]))
ax21 = fig.add_axes([0.55, 0.05, 0.35, 0.9])
plot2Ddata(
    receiver_locations[:, 0:2],
    dpred_plotting[:, -1],
    ax=ax21,
    ncontour=30,
    clim=(-v_max, v_max),
    contourOpts={"cmap": "bwr"},
)
ax21.set_title("dBz/dt at 0.001 s")

ax22 = fig.add_axes([0.92, 0.05, 0.02, 0.9])
norm2 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)
cbar2 = mpl.colorbar.ColorbarBase(
    ax22, norm=norm2, orientation="vertical", cmap=mpl.cm.bwr
)
cbar2.set_label("$T/s$", rotation=270, labelpad=15, size=12)

plt.show()

The variable dpred_plotting contains an array of the 121x31 of dB/dt. Plotting these results in clearly broken amplitudes and oscillatory behavior, and I’ve tried numerous ways to fix this.
I think I am not setting up the survey correctly (i.e. something to do with putting the source and receivers directly on the boundary). Does anyone have ideas?
Thanks much,

Julien

1 Like

Hi @jchaput, I took a quick look through. It looks like you are working with an interesting example! And I appreciated seeing the plotly code for 3D visualization.

I have a few suggestions that will hopefully be helpful here.

One important thing to note is that within SimPEG, we use fortran ordering, so for reshaping the data, you will want to use

dpred_plotting = np.reshape(dpred, (n_tx ** 2, n_times), order="F")

But an even more robust way is to use the Data class within SimPEG

from SimPEG import Data
data = Data(dobs=dpred, survey=survey)

and this allows you to then query by source / receiver pair

src = survey.source_list[0]
rx = src.receiver_list[0]
data[src, rx]

I coarsened the mesh to have it run quickly, and the values that are obtained are very large, so I am suspicious of the setup. On looking through it, there are a couple of things to be careful with…

Mesh

For an EM simulation, we need to include air cells in the simulation. It looks like the mesh designed here has a top at z=0m. This will cause issues as the fields need to have decayed at the boundaries, but here the source is located on that boundary. So this will cause issues.

Model & Mappings

It looks like the model is defined in terms of log-conductivity, but the mapping does not account for this (model = np.log(sig_background) * np.ones(mesh.nC) #background model). The mapping only considers injecting the active cells. An ExpMap will be necessary to convert the log-conductivity values to conductivity in units of S/m. Otherwise, the simulation is working with negative conductivity values, which will break. Currently, we don’t have safe-guards to make sure that the model is positive, but this would be a good thing for us to put a check in the code and at least give you a warning (warning or error for negative conductivity values in DC / EM simulations · Issue #1089 · simpeg/simpeg · GitHub).

Simulation & time-stepping

For the simulation, I would recommend starting with the step-off waveform. This is the simplest waveform and would let you test out the setup before introducing the complexity of a variable waveform

A couple of resources

A while back, I gave a seminar that was recorded on SimPEG, with a focus on EM. This goes over both the forward simulation and the inversion, so the first part if probably the most relevant to you. With this, I did a demo of a TDEM forward simulation and tried to walk through some of the important components such as mesh design. That notebook is here: time-domain-cyl-forward.ipynb - Curvenote. Note that some of the code might be behind the latest SimPEG, but the discussion is still relevant.

Finally, on the research connections-front, it working with highly conductive targets is definitely a topic of interest! We have done some work looking at steel-cased wells and recently published a short article in the leading edge with some discussion (doi, arxiv, notebooks). Anyways, happy to chat further on this!

Feel free to post again if you have further questions.

1 Like

Okay, I think I have a first version running (though I’m not convinced of the spatial patterns, the temporal decays and amplitudes look like they could be correct). Thanks for the Data object note; makes plotting a lot easier.
The forward model gives reasonable results if I don’t convert everything to log(conductivity) and don’t use an additional ExpMap, but if I try to implement the exp mapping, the results appear to be unstable, so I assume my implementation of the map is incorrect.
For instance this works:

air_value      = (1e-8)
sig_silver     = (1000)
sig_background = (2e-3)
model_map = (maps.InjectActiveCells(mesh, ind_active, air_value))

But this does not:

air_value      = np.log(1e-8)
sig_silver     = np.log(1000)
sig_background = np.log(2e-3)
model_map = (maps.InjectActiveCells(mesh, ind_active, air_value)* maps.ExpMap(nP=NN))

NN is the number of active cells (i.e. z<=0).
I have a few other questions and observations about the implementation of simpeg.
1- As you mentioned, we must define air cells even for a ground based simulation, but the conductivity model is only defined with respect to the active cells portion (I was getting all sort of array dimension problems for a while because I was defining the conductivity model with respect tot he whole mesh, and not the active cells). As such, when building arbitrary conductivity models, I found it useful to create a second mesh (meshG in the code) that is only used to reference cells to build the model.

2- The stepoff waveform currently yields believable dB/dt decay curves, but there are oddities. For instance, I currently have the time_off set to t0=-1e-5, and this is referenced in the simulation step as well. But, if I set the off time to t0=1e-5 (and match it in the simulation step), I get the "points outside of mesh error, even though my time_stepping only starts at 1e-5. Does this mean we can only time step the simulation for intervals AFTER the source has completely shut off?

3- On that note, the waveform I’m using in the field is a NanoTEM, so very rapid trapezoidal switch off. In terms of a trapezoid waveform, something like this:

waveform = tdem.sources.TrapezoidWaveform(
    ramp_on=np.r_[-7.81e-3, -7.8096e-3], ramp_off=np.r_[-0.0004e-3, 0.0], offTime=0.0
)

Given that it is trapezoidal but with rapid ramp on/off, should I simply sample the wave at its corners, i.e.:

waveform_times = np.array([-7.81e-3, -7.8096e-3,-0.0004e-3, 0.0])

Or should I densely and regularly sample it as in the other examples on the site? i.e.:

waveform_times = np.linspace(-7.81e-3,0,21)

4- My receivers also have an area (albeit smaller than the source loop) over which I also need to scale the readings. Is there a way to do that via options?

5- You mentioned that the field needs to have decayed at the boundaries (I assume there’s no PML or similar implemented here?). Do you have any intuition on how large the domain should be to ensure this?

Thanks! Here’s the code in one block for easy copy/paste (with the plotly commented out. Now also plots the source loop. I’ll post an easy to read one here with comments once it works.)

from SimPEG import Data
from discretize import TensorMesh
from discretize import TreeMesh
from discretize.utils import mkvc, refine_tree_xyz
import matplotlib.pyplot as plt
from SimPEG.utils import plot2Ddata, surface2ind_topo
from SimPEG import maps
import SimPEG.electromagnetics.time_domain as tdem
from SimPEG.utils import meshTensor
import numpy as np
import numpy.matlib as matlib
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
import plotly.offline as pyo
import plotly.graph_objs as go
import plotly.express as px
import plotly.io as pio
pio.renderers.default='browser'

%matplotlib inline
try:
    from pymatsolver import Pardiso as Solver
except ImportError:
    from SimPEG import SolverLU as Solver

save_file = False

'''
Create a mesh
'''

nc = 60  # number of core mesh cells in x, y and z
ncz = 28
dh = 20  # base cell width in x, y and z
dhz = 5
npad = 5  # number of padding cells: these are cells of varying size OUTSIDE the code cell region (usually less dense)
exp = 1.25  # expansion rate of padding cells

hx = [(dh, nc)] #reads like this: dilation left of zero, core cells, dilation positive from zero
hy = [(dh, nc)]
hz = [(dhz, 4, -exp), (dhz, ncz), (dhz, 4, exp)]
mesh = TensorMesh([hx, hy, hz], x0=[-600,-600,-150]) # creates the mesh
#mesh = TensorMesh([hx, hy, hz], x0="CCN") # creates the mesh


# The bottom southwest corner
x0 = mesh.x0 #indicates the southwest deepest cell

# The total number of cells
nC = mesh.nC

# An (nC, 3) array containing the cell-center locations
cc = mesh.gridCC

# A boolean array specifying which cells lie on the boundary
bInd = mesh.cellBoundaryInd

# The cell volumes
v = mesh.vol


ind_active = mesh.gridCC[:, 2] < 0.0
NN = np.count_nonzero(ind_active) #The number of active cells

# air_value      = np.log(1e-8)
# sig_silver     = np.log(1000)
# sig_background = np.log(2e-3)
# model_map = (maps.InjectActiveCells(mesh, ind_active, air_value)* maps.ExpMap(nP=NN))

air_value      = (1e-8)
sig_silver     = (1000)
sig_background = (2e-3)
model_map = (maps.InjectActiveCells(mesh, ind_active, air_value))#* maps.ExpMap(nP=NN))


#There a re subtleties with mapping. The forward model MUST have air cells defined, but the inverse model does not.
#The number of cells in the "mapping" must be the active cells
#NOTE: THE MODEL IS ONLY DEFINED ON ACTIVE CELLS, NOT THE WHOLE MESH


#Define a second mesh WITHOUT air cells?
meshG = TensorMesh([mesh.hx,mesh.hy, mesh.hz[mesh.vectorCCz <= 0]])

meshG.x0 = [-meshG.hx.sum() / 2.0, -meshG.hy.sum() / 2, -150]


'''
Define the conductivity model
'''

model = (sig_background) * np.ones(np.count_nonzero(ind_active)) #background model

#create line in 3D
xo,yo,zo = (0,0,-40) #passes through this point
x1,y1,z1 = (200,200,-65)
Lx,Ly,Lz = [np.linspace(xo,x1,100),np.linspace(yo,y1,100),np.linspace(zo,z1,100)]

#function to calculate the distance of every cell to a point
def dist(x,msh):
    dis = np.zeros(len(msh))
    p1 = np.array([meshG.gridCC[:,0], meshG.gridCC[:,1], meshG.gridCC[:,2]])
    p2 = np.transpose(matlib.repmat(np.transpose(x),len(p1[0]),1))
    dist = (((p2[0,:]-p1[0,:])**2)+((p2[1,:]-p1[1,:])**2)+((p2[2,:]-p1[2,:])**2))**(1/2)
    return dist

#Note: mesh.gridCC[x1,xyz] returns one element (x) and its 3d x,y,z position in the volume
thk = 20 #vein thickness

for i in range(len(Lx)):
    ds = dist(np.array([Lx[i],Ly[i],Lz[i]]),mesh)
    model[ds < thk] = (sig_silver)

'''
Define the receivers
'''

# Observation times for response (time channels)
n_times = 31
time_channels = np.logspace(-7, -3, n_times)

n_tx = 11
# Define receiver locations
xrx, yrx, zrx = np.meshgrid(
    np.linspace(-190, 190, n_tx), np.linspace(-190, 190, n_tx), [0]
)
receiver_locations = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)]
receivers_list = []

ntx = np.size(xrx)
# Each unique location defines a new transmitter
for ii in range(ntx):
    # Here we define receivers that measure the h-field in A/m
    dbzdt_receiver = tdem.receivers.PointMagneticFluxTimeDerivative(
        receiver_locations[ii, :], time_channels, orientation = "z"
    )
    receivers_list.append(dbzdt_receiver)
     # Make a list containing all receivers even if just one


'''
Define the source waveforms
'''
# waveform_times = np.linspace(-7.81e-3, 0, 51)

# # Set this up to emulate nanotem
# waveform = tdem.sources.TrapezoidWaveform(
#     ramp_on=np.r_[-7.81e-3, -7.8096e-3], ramp_off=np.r_[-0.0004e-3, 0.0], off_time=0.0
# )
# waveform_value = [waveform.eval(t) for t in waveform_times]

# waveform_times = np.linspace(-0.002, 0, 21)

# For each waveform type, you must define the necessary set of kwargs.
# For the trapezoidal waveform we define the ramp on interval, the
# ramp-off interval and the off-time.
# waveform = tdem.sources.TrapezoidWaveform(
#     ramp_on=np.r_[-0.007, -0.006], ramp_off=np.r_[-0.001, 0.0], off_time=0.0
# )

waveform = tdem.sources.StepOffWaveform(off_time=-1e-5)


# waveform_value = [waveform.eval(t) for t in waveform_times]



'''
plot the survey
'''

#define source loop on a circle (densely sampled loop)
# angs = np.arange(0.01,2*np.pi,0.05)
# RD = 28
# xs = RD*np.cos(angs)
# ys = RD*np.sin(angs)
# zs = np.zeros((1,len(xs)))
# SL = np.c_[mkvc(xs), mkvc(ys), mkvc(zs)]


# fig =go.Figure(data=go.Isosurface(
#     x= meshG.gridCC[:,0],
#     y= meshG.gridCC[:,1],
#     z= meshG.gridCC[:,2],
#     value=model,
#     isomin=500,
#     isomax=1500,
# ))

# fig.add_scatter3d(    
#     x=receiver_locations[:,0],
#     y=receiver_locations[:,1],
#     z=receiver_locations[:,2],
#     mode='markers',
# )
# fig.add_scatter3d(    
#     x=SL[:,0],
#     y=SL[:,1],
#     z=SL[:,2],
#     mode='markers',

# )


# fig.show()



'''
Define the source
'''

src_ramp_on = tdem.Src.CircularLoop(
    receivers_list,
    location=np.r_[0.0, 0.0, 0.0],
    orientation = "z",
    current = 2.5,
    radius=28,
    waveform=waveform,
)
src_list_ramp_on = [src_ramp_on]

survey_ramp_on = tdem.Survey(src_list_ramp_on)

survey = tdem.Survey(src_list_ramp_on)

'''
Simulation
'''

# dt and how many dt will result (so here we have 40 time steps in total)
time_steps = [(1e-5, 20), (1e-5, 10), (1e-4, 30), (1e-3, 5)]

#check total time
dt = meshTensor(time_steps)
tmax = np.hstack([[0], dt]).cumsum()

#Note: waveform on time was defined as -0.002 to 0
simulation = tdem.simulation.Simulation3DMagneticFluxDensity(
    mesh, survey=survey, sigmaMap=model_map,verbose=True, solver=Solver, t0=-1e-5
)

# Set the time-stepping for the simulation
simulation.time_steps = time_steps

dpred = simulation.dpred(model)

'''
plot things
'''

#alternate plotting/referencing code : plot individual dB/dt decays
data = Data(dobs=dpred, survey=survey)
src = survey.source_list[0]
for i in range(5):
    rx = src.receiver_list[i]
    plt.plot(data[src, rx])

# Data were organized by location, then by time channel
dpred_plotting = np.reshape(dpred, (n_tx ** 2, n_times), order="F")

# Plot
fig = plt.figure(figsize=(10, 4))

# dB/dt at early time
v_max = np.max(np.abs(dpred_plotting[:, 0]))
ax11 = fig.add_axes([0.05, 0.05, 0.35, 0.9])
plot2Ddata(
    receiver_locations[:, 0:2],
    dpred_plotting[:, 0],
    ax=ax11,
    ncontour=30,
    clim=(-v_max, v_max),
    contourOpts={"cmap": "bwr"},
)
ax11.set_title("dBz/dt at 0.0001 s")

ax12 = fig.add_axes([0.42, 0.05, 0.02, 0.9])
norm1 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)
cbar1 = mpl.colorbar.ColorbarBase(
    ax12, norm=norm1, orientation="vertical", cmap=mpl.cm.bwr
)
cbar1.set_label("$T/s$", rotation=270, labelpad=15, size=12)


# dB/dt at later times
v_max = np.max(np.abs(dpred_plotting[:, -1]))
ax21 = fig.add_axes([0.55, 0.05, 0.35, 0.9])
plot2Ddata(
    receiver_locations[:, 0:2],
    dpred_plotting[:, -1],
    ax=ax21,
    ncontour=30,
    clim=(-v_max, v_max),
    contourOpts={"cmap": "bwr"},
)
ax21.set_title("dBz/dt at 0.001 s")

ax22 = fig.add_axes([0.92, 0.05, 0.02, 0.9])
norm2 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)
cbar2 = mpl.colorbar.ColorbarBase(
    ax22, norm=norm2, orientation="vertical", cmap=mpl.cm.bwr
)
cbar2.set_label("$T/s$", rotation=270, labelpad=15, size=12)

plt.show()