If you a have a series of 2D rasters with the same shape, you can stack them to a 3D array:
import numpy as np
Z = np.stack((Z1, Z2, Z3, ...))
Then use min, max or similar along the first axis (i.e.
Zmin = Z.min(0)
Zmax = Z.max(0)
Which should result in a 2D array the same shape as one of the input arrays.
For very large files, there are a few ideas. If your 2D rasters are larger than your available RAM, you will need to do multiple passes with window subsets (nXOff, nYOff, nXSize, nYSize) on each raster, which you can do with GDAL. Each window stack would find the minimum, and write the result.
However, if you can hold a few 2D rasters in your RAM at once, you could read them one-by-one and accumulate the result:
shape = (10000, 10000)
Zmin = np.inf * np.ones(shape, 'f') # Float32
Zmax = -np.inf * np.ones(shape, 'f')
for file in file_list:
ds = gdal.Open(file)
ar = ds.ReadAsArray()
ds = None # close
assert ar.shape == shape
sel = ar < Zmin
Zmin[sel] = ar[sel]
sel = ar > Zmax
Zmax[sel] = ar[sel]