docs: def discharge(self):
          """Discharge of the element in each layer
          """
          
          rv = np.zeros(self.aq[0].naq)
          Qls = self.parameters[:, 0] * self.dischargeinf()
          Qls.shape = (self.nls, self.nlayers, self.order + 1)
          Qls = np.sum(Qls, 2)
          for i, q in enumerate(Qls):
              rv[self.layers[i]] += q
              #rv[self.layers] = np.sum(Qls.reshape(self.nls * (self.order + 1), self.nlayers), 0)
          return rv