There are several approaches to answer this question. One is simple, but not that accurate. This approach would be to to discretize points along the perimeter and calculate a per point boundary shear stress. The equation would be the same as the depth averaged, except you would solve for each individual depth h1, ..., hn.
One thing to note is that the bed load will mostly be concentrated at the lowest point in the channel.
The easiest accurate method to implement is the WTA method, which uses the law of the wall to calculate boundary shear stress. The papers outlining the method are Wobus et al. (2006, 2008) titled "Self-formed bedrock channels", and "Modeling the evolution of channel shape: Balancing computational efficiency with hydraulic fidelity" . I have a python code that implements this, but in closed conduit conditions. It can be modified to work in the open channel case.
- Wobus et al., 2006: Self-formed bedrock channels. doi: 10.1029/2006GL027182
- Wobus et al., 2008: Modeling the evolution of channel shape: Balancing computational efficiency with hydraulic fidelity. doi: 10.1029/2007JF000914