R code for a constrained randomization in a three-dimensional array?

Starting with this:

set.seed(42)
ary  apply(1, mean)
        A         B         C 
10.857143 10.857143  6.285714 

However, I think you may be under-constraining here: while this meets your "average" constraints, it may not meet reasonable distributions:

### Each pipe ...
apply(ary, 1:2, table) |> apply(1, range)
      A B C
[1,]  1 1 1
[2,] 10 9 6
### Each depth level ...
apply(ary, 3, table) |> apply(1, range)
      A  B C
[1,]  9  6 3
[2,] 15 14 9

This means that while you do average around 5.4 for each of A and B, they both range from 1-9 within any given pipe.

One can work around this by randomly iterating. Let's say you'd be willing to accept ranges of 3-8 for A and B; similarly, enforcing ranges of (say) 8-16 for the depth-level tables.

Side note: the apply code worked above, but occasionally when one or more of the treatments are missing in a particular pipe, then it won't work. I'll create a table2 that is more immune to this problem.

set.seed(43)
ary 

This is excellent! I think I was unclear in the specification of the constraints, though: ideally the number of samples per treatment should be either 5 or 6 per X,Y coordinate for the treatments with more reps, and 3 or 4 for the treatments with fewer. That is reflected in the code I edited in to the answer.